diff --git a/kwant/linalg/decomp_ev.py b/kwant/linalg/decomp_ev.py
index bd92ce9335d2ba31edc4535f0a8a12a5828e2943..8e7b95d68b0eae6701f3fe38a624653a5aa67f68 100644
--- a/kwant/linalg/decomp_ev.py
+++ b/kwant/linalg/decomp_ev.py
@@ -50,18 +50,5 @@ def gen_eig(a, b, left=False, right=True, overwrite_ab=False):
The right eigenvector corresponding to the eigenvalue
``alpha[i]/beta[i]`` is the column ``vr[:,i]``.
"""
-
- ltype, a, b = lapack.prepare_for_lapack(overwrite_ab, a, b)
-
- if a.ndim != 2 or b.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 b.shape[0] != a.shape[0] or b.shape[1] != a.shape[1]:
- raise ValueError("gen_eig requires a and be to have the same shape")
-
- ggev = getattr(lapack, ltype + "ggev")
-
- return ggev(a, b, left, right)
+ a, b = lapack.prepare_for_lapack(overwrite_ab, a, b)
+ return lapack.ggev(a, b, left, right)
diff --git a/kwant/linalg/decomp_lu.py b/kwant/linalg/decomp_lu.py
index 79accc254feaa53206a9ea2f4a9d977a0d904138..639b547bfdb082eb0e45df3180bb31ceeb67937f 100644
--- a/kwant/linalg/decomp_lu.py
+++ b/kwant/linalg/decomp_lu.py
@@ -45,20 +45,8 @@ def lu_factor(a, overwrite_a=False):
singular : boolean
Whether the matrix a is singular (up to machine precision)
"""
-
- ltype, a = lapack.prepare_for_lapack(overwrite_a, a)
-
- if a.ndim != 2:
- raise ValueError("lu_factor expects a matrix")
-
- if ltype == 'd':
- return lapack.dgetrf(a)
- elif ltype == 'z':
- return lapack.zgetrf(a)
- elif ltype == 's':
- return lapack.sgetrf(a)
- else:
- return lapack.cgetrf(a)
+ a = lapack.prepare_for_lapack(overwrite_a, a)
+ return lapack.getrf(a)
def lu_solve(matrix_factorization, b):
@@ -83,23 +71,9 @@ def lu_solve(matrix_factorization, b):
"a singular matrix. Result of solve step "
"are probably unreliable")
- ltype, lu, b = lapack.prepare_for_lapack(False, lu, b)
+ lu, b = lapack.prepare_for_lapack(False, lu, b)
ipiv = np.ascontiguousarray(np.asanyarray(ipiv), dtype=lapack.int_dtype)
-
- if b.ndim > 2:
- raise ValueError("lu_solve: b must be a vector or matrix")
-
- if lu.shape[0] != b.shape[0]:
- raise ValueError("lu_solve: incompatible dimensions of b")
-
- if ltype == 'd':
- return lapack.dgetrs(lu, ipiv, b)
- elif ltype == 'z':
- return lapack.zgetrs(lu, ipiv, b)
- elif ltype == 's':
- return lapack.sgetrs(lu, ipiv, b)
- else:
- return lapack.cgetrs(lu, ipiv, b)
+ return lapack.getrs(lu, ipiv, b)
def rcond_from_lu(matrix_factorization, norm_a, norm="1"):
@@ -127,17 +101,6 @@ def rcond_from_lu(matrix_factorization, norm_a, norm="1"):
norm specified in norm
"""
(lu, ipiv, singular) = matrix_factorization
- if not norm in ("1", "I"):
- raise ValueError("norm in rcond_from_lu must be either '1' or 'I'")
norm = norm.encode('utf8') # lapack expects bytes
-
- ltype, lu = lapack.prepare_for_lapack(False, lu)
-
- if ltype == 'd':
- return lapack.dgecon(lu, norm_a, norm)
- elif ltype == 'z':
- return lapack.zgecon(lu, norm_a, norm)
- elif ltype == 's':
- return lapack.sgecon(lu, norm_a, norm)
- else:
- return lapack.cgecon(lu, norm_a, norm)
+ lu = lapack.prepare_for_lapack(False, lu)
+ return lapack.gecon(lu, norm_a, norm)
diff --git a/kwant/linalg/decomp_schur.py b/kwant/linalg/decomp_schur.py
index 4accc17ff0f1643c9532ada85a29dfa08a2b67fa..d65432132efd3084e26a25643dbc36b82f48b39c 100644
--- a/kwant/linalg/decomp_schur.py
+++ b/kwant/linalg/decomp_schur.py
@@ -62,18 +62,8 @@ def schur(a, calc_q=True, calc_ev=True, overwrite_a=False):
LinAlgError
If the underlying QR iteration fails to converge.
"""
-
- ltype, a = lapack.prepare_for_lapack(overwrite_a, a)
-
- if a.ndim != 2:
- raise ValueError("Expect matrix as input")
-
- if a.shape[0] != a.shape[1]:
- raise ValueError("Expect square matrix")
-
- gees = getattr(lapack, ltype + "gees")
-
- return gees(a, calc_q, calc_ev)
+ a = lapack.prepare_for_lapack(overwrite_a, a)
+ return lapack.gees(a, calc_q, calc_ev)
def convert_r2c_schur(t, q):
@@ -192,9 +182,7 @@ def order_schur(select, t, q, calc_ev=True, overwrite_tq=False):
``calc_ev == True``
"""
- ltype, t, q = lapack.prepare_for_lapack(overwrite_tq, t, q)
-
- trsen = getattr(lapack, ltype + "trsen")
+ t, q = lapack.prepare_for_lapack(overwrite_tq, t, q)
# Figure out if select is a function or array.
isfun = isarray = True
@@ -223,7 +211,7 @@ def order_schur(select, t, q, calc_ev=True, overwrite_tq=False):
t, q = convert_r2c_schur(t, q)
return order_schur(select, t, q, calc_ev, True)
- return trsen(select, t, q, calc_ev)
+ return lapack.trsen(select, t, q, calc_ev)
def evecs_from_schur(t, q, select=None, left=False, right=True,
@@ -267,13 +255,7 @@ def evecs_from_schur(t, q, select=None, left=False, right=True,
``right == True``.
"""
- ltype, t, q = lapack.prepare_for_lapack(overwrite_tq, t, q)
-
- if (t.shape[0] != t.shape[1] or q.shape[0] != q.shape[1]
- or t.shape[0] != q.shape[0]):
- raise ValueError("Invalid Schur decomposition as input")
-
- trevc = getattr(lapack, ltype + "trevc")
+ t, q = lapack.prepare_for_lapack(overwrite_tq, t, q)
# check if select is a function or an array
if select is not None:
@@ -300,7 +282,7 @@ def evecs_from_schur(t, q, select=None, left=False, right=True,
else:
selectarr = None
- return trevc(t, q, selectarr, left, right)
+ return lapack.trevc(t, q, selectarr, left, right)
def gen_schur(a, b, calc_q=True, calc_z=True, calc_ev=True,
@@ -364,21 +346,8 @@ def gen_schur(a, b, calc_q=True, calc_z=True, calc_ev=True,
LinAlError
If the underlying QZ iteration fails to converge.
"""
-
- ltype, a, b = lapack.prepare_for_lapack(overwrite_ab, a, b)
-
- if a.ndim != 2 or b.ndim != 2:
- raise ValueError("Expect matrices as input")
-
- if a.shape[0] != a.shape[1]:
- raise ValueError("Expect square matrix a")
-
- if a.shape[0] != b.shape[0] or a.shape[0] != b.shape[1]:
- raise ValueError("Shape of b is incompatible to matrix a")
-
- gges = getattr(lapack, ltype + "gges")
-
- return gges(a, b, calc_q, calc_z, calc_ev)
+ a, b = lapack.prepare_for_lapack(overwrite_ab, a, b)
+ return lapack.gges(a, b, calc_q, calc_z, calc_ev)
def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
@@ -439,22 +408,8 @@ def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
LinAlError
If the problem is too ill-conditioned.
"""
- ltype, s, t, q, z = lapack.prepare_for_lapack(overwrite_stqz, s, t, q, z)
+ s, t, q, z = lapack.prepare_for_lapack(overwrite_stqz, s, t, q, z)
- if (s.ndim != 2 or t.ndim != 2 or
- (q is not None and q.ndim != 2) or
- (z is not None and z.ndim != 2)):
- raise ValueError("Expect matrices as input")
-
- if ((s.shape[0] != s.shape[1] or t.shape[0] != t.shape[1] or
- s.shape[0] != t.shape[0]) or
- (q is not None and (q.shape[0] != q.shape[1] or
- s.shape[0] != q.shape[0])) or
- (z is not None and (z.shape[0] != z.shape[1] or
- s.shape[0] != z.shape[0]))):
- raise ValueError("Invalid Schur decomposition as input")
-
- tgsen = getattr(lapack, ltype + "tgsen")
# Figure out if select is a function or array.
isfun = isarray = True
@@ -492,7 +447,7 @@ def order_gen_schur(select, s, t, q=None, z=None, calc_ev=True,
return order_gen_schur(select, s, t, q, z, calc_ev, True)
- return tgsen(select, s, t, q, z, calc_ev)
+ return lapack.tgsen(select, s, t, q, z, calc_ev)
def convert_r2c_gen_schur(s, t, q=None, z=None):
@@ -536,7 +491,7 @@ def convert_r2c_gen_schur(s, t, q=None, z=None):
If it fails to convert a 2x2 block into complex form (unlikely).
"""
- ltype, s, t, q, z = lapack.prepare_for_lapack(True, s, t, q, z)
+ s, t, q, z = lapack.prepare_for_lapack(True, s, t, q, z)
# Note: overwrite=True does not mean much here, the arrays are all copied
if (s.ndim != 2 or t.ndim != 2 or
@@ -656,20 +611,7 @@ def evecs_from_gen_schur(s, t, q=None, z=None, select=None,
"""
- ltype, s, t, q, z = lapack.prepare_for_lapack(overwrite_qz, s, t, q, z)
-
- if (s.ndim != 2 or t.ndim != 2 or
- (q is not None and q.ndim != 2) or
- (z is not None and z.ndim != 2)):
- raise ValueError("Expect matrices as input")
-
- if ((s.shape[0] != s.shape[1] or t.shape[0] != t.shape[1] or
- s.shape[0] != t.shape[0]) or
- (q is not None and (q.shape[0] != q.shape[1] or
- s.shape[0] != q.shape[0])) or
- (z is not None and (z.shape[0] != z.shape[1] or
- s.shape[0] != z.shape[0]))):
- raise ValueError("Invalid Schur decomposition as input")
+ s, t, q, z = lapack.prepare_for_lapack(overwrite_qz, s, t, q, z)
if left and q is None:
raise ValueError("Matrix q must be provided for left eigenvectors")
@@ -677,8 +619,6 @@ def evecs_from_gen_schur(s, t, q=None, z=None, select=None,
if right and z is None:
raise ValueError("Matrix z must be provided for right eigenvectors")
- tgevc = getattr(lapack, ltype + "tgevc")
-
# Check if select is a function or an array.
if select is not None:
isfun = isarray = True
@@ -704,4 +644,4 @@ def evecs_from_gen_schur(s, t, q=None, z=None, select=None,
else:
selectarr = None
- return tgevc(s, t, q, z, selectarr, left, right)
+ return lapack.tgevc(s, t, q, z, selectarr, left, right)
diff --git a/kwant/linalg/f_lapack.pxd b/kwant/linalg/f_lapack.pxd
deleted file mode 100644
index 1a0304af871c6bb75994eb0668389e0d6ad2b54d..0000000000000000000000000000000000000000
--- a/kwant/linalg/f_lapack.pxd
+++ /dev/null
@@ -1,217 +0,0 @@
-# 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
-# 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.
-
-ctypedef int l_int
-ctypedef int l_logical
-
-cdef extern:
- void sgetrf_(l_int *, l_int *, float *, l_int *, l_int *, l_int *)
- void dgetrf_(l_int *, l_int *, double *, l_int *, l_int *, l_int *)
- void cgetrf_(l_int *, l_int *, float complex *, l_int *, l_int *,
- l_int *)
- void zgetrf_(l_int *, l_int *, double complex *, l_int *, l_int *,
- l_int *)
-
- void sgetrs_(char *, l_int *, l_int *, float *, l_int *, l_int *,
- float *, l_int *, l_int *)
- void dgetrs_(char *, l_int *, l_int *, double *, l_int *, l_int *,
- double *, l_int *, l_int *)
- void cgetrs_(char *, l_int *, l_int *, float complex *, l_int *,
- l_int *, float complex *, l_int *, l_int *)
- void zgetrs_(char *, l_int *, l_int *, double complex *, l_int *,
- l_int *, double complex *, l_int *, l_int *)
-
- void sgecon_(char *, l_int *, float *, l_int *, float *, float *,
- float *, l_int *, l_int *)
- void dgecon_(char *, l_int *, double *, l_int *, double *, double *,
- double *, l_int *, l_int *)
- void cgecon_(char *, l_int *, float complex *, l_int *, float *,
- float *, float complex *, float *, l_int *)
- void zgecon_(char *, l_int *, double complex *, l_int *, double *,
- double *, double complex *, double *, l_int *)
-
- void sggev_(char *, char *, l_int *, float *, l_int *, float *, l_int *,
- float *, float *, float *, float *, l_int *, float *, l_int *,
- float *, l_int *, l_int *)
- void dggev_(char *, char *, l_int *, double *, l_int *, double *, l_int *,
- double *, double *, double *, double *, l_int *,
- double *, l_int *, double *, l_int *, l_int *)
- void cggev_(char *, char *, l_int *, float complex *, l_int *,
- float complex *, l_int *, float complex *, float complex *,
- float complex *, l_int *, float complex *, l_int *,
- float complex *, l_int *, float *, l_int *)
- void zggev_(char *, char *, l_int *, double complex *, l_int *,
- double complex *, l_int *, double complex *,
- double complex *, double complex *, l_int *,
- double complex *, l_int *, double complex *, l_int *,
- double *, l_int *)
-
- void sgees_(char *, char *, l_logical (*)(float *, float *),
- l_int *, float *, l_int *, l_int *,
- float *, float *, float *, l_int *,
- float *, l_int *, l_logical *, l_int *)
- void dgees_(char *, char *, l_logical (*)(double *, double *),
- l_int *, double *, l_int *, l_int *,
- double *, double *, double *, l_int *,
- double *, l_int *, l_logical *, l_int *)
- void cgees_(char *, char *,
- l_logical (*)(float complex *),
- l_int *, float complex *,
- l_int *, l_int *, float complex *,
- float complex *, l_int *,
- float complex *, l_int *, float *,
- l_logical *, l_int *)
- void zgees_(char *, char *,
- l_logical (*)(double complex *),
- l_int *, double complex *,
- l_int *, l_int *, double complex *,
- double complex *, l_int *,
- double complex *, l_int *,
- double *, l_logical *, l_int *)
-
- void strsen_(char *, char *, l_logical *, l_int *,
- float *, l_int *, float *,
- l_int *, float *, float *, l_int *,
- float *, float *, float *, l_int *,
- l_int *, l_int *, l_int *)
- void dtrsen_(char *, char *, l_logical *,
- l_int *, double *, l_int *,
- double *, l_int *, double *, double *,
- l_int *, double *, double *, double *,
- l_int *, l_int *, l_int *, l_int *)
- void ctrsen_(char *, char *, l_logical *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, float complex *, l_int *,
- float *, float *, float complex *,
- l_int *, l_int *)
- void ztrsen_(char *, char *, l_logical *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, double complex *, l_int *,
- double *, double *, double complex *,
- l_int *, l_int *)
-
- void strevc_(char *, char *, l_logical *,
- l_int *, float *, l_int *,
- float *, l_int *, float *, l_int *,
- l_int *, l_int *, float *, l_int *)
- void dtrevc_(char *, char *, l_logical *,
- l_int *, double *, l_int *,
- double *, l_int *, double *,
- l_int *, l_int *, l_int *, double *,
- l_int *)
- void ctrevc_(char *, char *, l_logical *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, l_int *, l_int *,
- float complex *, float *, l_int *)
- void ztrevc_(char *, char *, l_logical *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, l_int *, l_int *,
- double complex *, double *, l_int *)
-
- void sgges_(char *, char *, char *,
- l_logical (*)(float *, float *, float *),
- l_int *, float *, l_int *, float *,
- l_int *, l_int *, float *, float *,
- float *, float *, l_int *, float *,
- l_int *, float *, l_int *, l_logical *,
- l_int *)
- void dgges_(char *, char *, char *,
- l_logical (*)(double *, double *, double *),
- l_int *, double *, l_int *, double *,
- l_int *, l_int *, double *, double *,
- double *, double *, l_int *, double *,
- l_int *, double *, l_int *,
- l_logical *, l_int *)
- void cgges_(char *, char *, char *,
- l_logical (*)(float complex *, float complex *),
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, l_int *, float complex *,
- float complex *, float complex *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, float *, l_logical *, l_int *)
- void zgges_(char *, char *, char *,
- l_logical (*)(double complex *, double complex *),
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, l_int *, double complex *,
- double complex *,
- double complex *, l_int *,
- double complex *, l_int *,
- double complex *, l_int *,
- double *, l_logical *, l_int *)
-
- void stgsen_(l_int *, l_logical *,
- l_logical *, l_logical *,
- l_int *, float *, l_int *, float *,
- l_int *, float *, float *, float *,
- float *, l_int *, float *, l_int *,
- l_int *, float *, float *, float *, float *,
- l_int *, l_int *, l_int *, l_int *)
- void dtgsen_(l_int *, l_logical *,
- l_logical *, l_logical *,
- l_int *, double *, l_int *,
- double *, l_int *, double *, double *,
- double *, double *, l_int *, double *,
- l_int *, l_int *, double *, double *,
- double *, double *, l_int *, l_int *,
- l_int *, l_int *)
- void ctgsen_(l_int *, l_logical *,
- l_logical *, l_logical *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, float complex *,
- float complex *,
- float complex *, l_int *,
- float complex *, l_int *, l_int *,
- float *, float *, float *,
- float complex *, l_int *, l_int *,
- l_int *, l_int *)
- void ztgsen_(l_int *, l_logical *,
- l_logical *, l_logical *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, double complex *,
- double complex *,
- double complex *, l_int *,
- double complex *, l_int *, l_int *,
- double *, double *, double *,
- double complex *, l_int *, l_int *,
- l_int *, l_int *)
-
- void stgevc_(char *, char *, l_logical *,
- l_int *, float *, l_int *,
- float *, l_int *, float *,
- l_int *, float *, l_int *,
- l_int *, l_int *, float *, l_int *)
- void dtgevc_(char *, char *, l_logical *,
- l_int *, double *, l_int *,
- double *, l_int *, double *,
- l_int *, double *, l_int *,
- l_int *, l_int *, double *, l_int *)
- void ctgevc_(char *, char *, l_logical *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, float complex *,
- l_int *, l_int *, l_int *,
- float complex *, float *, l_int *)
- void ztgevc_(char *, char *, l_logical *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, double complex *,
- l_int *, l_int *, l_int *,
- double complex *, double *, l_int *)
diff --git a/kwant/linalg/f_lapack.py b/kwant/linalg/f_lapack.py
deleted file mode 100644
index 39b702b5394dafb773744dfa1d0faa05c17a419c..0000000000000000000000000000000000000000
--- a/kwant/linalg/f_lapack.py
+++ /dev/null
@@ -1,14 +0,0 @@
-# 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
-# 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.
-
-__all__ = ['l_int_dtype', 'l_logical_dtype']
-
-import numpy as np
-
-l_int_dtype = np.int32
-l_logical_dtype = np.int32
diff --git a/kwant/linalg/lapack.pyx b/kwant/linalg/lapack.pyx
index d6699a0e11ea912506d417d1e1ab59270c3270a3..e173bc0ce28e860f93ee58d36c725d020fd95c72 100644
--- a/kwant/linalg/lapack.pyx
+++ b/kwant/linalg/lapack.pyx
@@ -1,4 +1,4 @@
-# Copyright 2011-2013 Kwant authors.
+# Copyright 2011-2017 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
@@ -8,29 +8,53 @@
"""Low-level access to LAPACK functions. """
-__all__ = ['sgetrf', 'dgetrf', 'cgetrf', 'zgetrf',
- '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',
+__all__ = ['getrf',
+ 'getrs',
+ 'gecon',
+ 'ggev',
+ 'gees',
+ 'trsen',
+ 'trevc',
+ 'gges',
+ 'tgsen',
+ 'tgevc',
'prepare_for_lapack']
import numpy as np
cimport numpy as np
-# Strangely absolute imports from kwant.linalg.f_lapack don't work here. So
-# continue using traditional implicit relative imports.
-from . import f_lapack
-from . cimport f_lapack
-from .f_lapack cimport l_int, l_logical
+cimport scipy.linalg.cython_lapack as lapack
-int_dtype = f_lapack.l_int_dtype
-logical_dtype = f_lapack.l_logical_dtype
+ctypedef int l_int
+ctypedef bint l_logical
+
+int_dtype = np.int32
+logical_dtype = np.int32
+
+ctypedef float complex float_complex
+ctypedef double complex double_complex
+
+ctypedef fused scalar:
+ float
+ double
+ float complex
+ double complex
+
+ctypedef fused single_precision:
+ float
+ float complex
+
+ctypedef fused double_precision:
+ double
+ double complex
+
+ctypedef fused cmplx:
+ float complex
+ double complex
+
+ctypedef fused floating:
+ float
+ double
# exceptions
@@ -51,284 +75,161 @@ def assert_fortran_mat(*mats):
raise ValueError("Input matrix must be Fortran contiguous")
-# Wrappers for xGETRF
-def sgetrf(np.ndarray[np.float32_t, ndim=2] A):
- cdef l_int M, N, info
- cdef np.ndarray[l_int, ndim=1] ipiv
-
- assert_fortran_mat(A)
-
- M = A.shape[0]
- N = A.shape[1]
- ipiv = np.empty(min(M,N), dtype = f_lapack.l_int_dtype)
-
- f_lapack.sgetrf_(&M, &N, <float *>A.data, &M,
- <l_int *>ipiv.data, &info)
-
- assert info >= 0, "Argument error in sgetrf"
-
- return (A, ipiv, info > 0 or M != N)
-
-def dgetrf(np.ndarray[np.float64_t, ndim=2] A):
- cdef l_int M, N, info
- cdef np.ndarray[l_int, ndim=1] ipiv
-
- assert_fortran_mat(A)
-
- M = A.shape[0]
- N = A.shape[1]
- ipiv = np.empty(min(M,N), dtype = f_lapack.l_int_dtype)
-
- f_lapack.dgetrf_(&M, &N, <double *>A.data, &M,
- <l_int *>ipiv.data, &info)
-
- assert info >= 0, "Argument error in dgetrf"
-
- return (A, ipiv, info > 0 or M != N)
-
-def cgetrf(np.ndarray[np.complex64_t, ndim=2] A):
- cdef l_int M, N, info
- cdef np.ndarray[l_int, ndim=1] ipiv
-
- assert_fortran_mat(A)
-
- M = A.shape[0]
- N = A.shape[1]
- ipiv = np.empty(min(M,N), dtype = f_lapack.l_int_dtype)
+cdef np.ndarray maybe_complex(scalar selector,
+ np.ndarray real, np.ndarray imag):
+ cdef np.ndarray r
+ r = real
+ if scalar in floating:
+ if imag.nonzero()[0].size:
+ r = real + 1j * imag
+ return r
- f_lapack.cgetrf_(&M, &N, <float complex *>A.data, &M,
- <l_int *>ipiv.data, &info)
- assert info >= 0, "Argument error in cgetrf"
+cdef l_int lwork_from_qwork(scalar qwork):
+ if scalar in floating:
+ return <l_int>qwork
+ else:
+ return <l_int>qwork.real
- return (A, ipiv, info > 0 or M != N)
-def zgetrf(np.ndarray[np.complex128_t, ndim=2] A):
+def getrf(np.ndarray[scalar, ndim=2] A):
cdef l_int M, N, info
- cdef np.ndarray[l_int, ndim=1] ipiv
+ cdef np.ndarray[l_int] ipiv
assert_fortran_mat(A)
M = A.shape[0]
N = A.shape[1]
- ipiv = np.empty(min(M,N), dtype = f_lapack.l_int_dtype)
-
- f_lapack.zgetrf_(&M, &N, <double complex *>A.data, &M,
- <l_int *>ipiv.data, &info)
-
- assert info >= 0, "Argument error in zgetrf"
+ ipiv = np.empty(min(M,N), dtype = int_dtype)
+
+ if scalar is float:
+ lapack.sgetrf(&M, &N, <float *>A.data, &M,
+ <l_int *>ipiv.data, &info)
+ elif scalar is double:
+ lapack.dgetrf(&M, &N, <double *>A.data, &M,
+ <l_int *>ipiv.data, &info)
+ elif scalar is float_complex:
+ lapack.cgetrf(&M, &N, <float complex *>A.data, &M,
+ <l_int *>ipiv.data, &info)
+ elif scalar is double_complex:
+ lapack.zgetrf(&M, &N, <double complex *>A.data, &M,
+ <l_int *>ipiv.data, &info)
+
+ assert info >= 0, "Argument error in getrf"
return (A, ipiv, info > 0 or M != N)
-# Wrappers for xGETRS
-def sgetrs(np.ndarray[np.float32_t, ndim=2] LU,
- np.ndarray[l_int, ndim=1] IPIV, B):
+def getrs(np.ndarray[scalar, ndim=2] LU, np.ndarray[l_int] IPIV,
+ np.ndarray B):
cdef l_int N, NRHS, info
- cdef np.ndarray b
assert_fortran_mat(LU)
- # again: workaround for 1x1-Fortran bug in NumPy < v2.0
- if (not isinstance(B, np.ndarray) or
- (B.ndim == 2 and (B.shape[0] > 1 or B.shape[1] > 1) and
- not B.flags["F_CONTIGUOUS"])):
- raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
-
- b = B
- N = LU.shape[0]
- if b.ndim == 1:
- NRHS = 1
- elif b.ndim == 2:
- NRHS = B.shape[1]
- else:
- raise ValueError("In sgetrs: B must be a vector or matrix")
-
- f_lapack.sgetrs_("N", &N, &NRHS, <float *>LU.data, &N,
- <l_int *>IPIV.data, <float *>b.data, &N,
- &info)
-
- assert info == 0, "Argument error in sgetrs"
-
- return b
+ # Consistency checks for LU and B
-def dgetrs(np.ndarray[np.float64_t, ndim=2] LU,
- np.ndarray[l_int, ndim=1] IPIV, B):
- cdef l_int N, NRHS, info
- cdef np.ndarray b
-
- assert_fortran_mat(LU)
+ if B.descr.type_num != LU.descr.type_num:
+ raise TypeError('B must have same dtype as LU')
- # again: workaround for 1x1-Fortran bug in NumPy < v2.0
- if (not isinstance(B, np.ndarray) or
- (B.ndim == 2 and (B.shape[0] > 1 or B.shape[1] > 1) and
+ # Workaround for 1x1-Fortran bug in NumPy < v2.0
+ if ((B.ndim == 2 and (B.shape[0] > 1 or B.shape[1] > 1) and
not B.flags["F_CONTIGUOUS"])):
- raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
-
- b = B
- N = LU.shape[0]
- if b.ndim == 1:
- NRHS = 1
- elif b.ndim == 2:
- NRHS = b.shape[1]
- else:
- raise ValueError("In dgetrs: B must be a vector or matrix")
-
- f_lapack.dgetrs_("N", &N, &NRHS, <double *>LU.data, &N,
- <l_int *>IPIV.data, <double *>b.data, &N,
- &info)
+ raise ValueError("B must be Fortran ordered")
- assert info == 0, "Argument error in dgetrs"
+ if B.ndim > 2:
+ raise ValueError("B must be a vector or matrix")
- return b
+ if LU.shape[0] != B.shape[0]:
+ raise ValueError('LU and B have incompatible shapes')
-def cgetrs(np.ndarray[np.complex64_t, ndim=2] LU,
- np.ndarray[l_int, ndim=1] IPIV, B):
- cdef l_int N, NRHS, info
- cdef np.ndarray b
-
- assert_fortran_mat(LU)
-
- # again: workaround for 1x1-Fortran bug in NumPy < v2.0
- if (not isinstance(B, np.ndarray) or
- (B.ndim == 2 and (B.shape[0] > 1 or B.shape[1] > 1) and
- not B.flags["F_CONTIGUOUS"])):
- raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
-
- b = B
N = LU.shape[0]
- if b.ndim == 1:
- NRHS = 1
- elif b.ndim == 2:
- NRHS = b.shape[1]
- else:
- raise ValueError("In cgetrs: B must be a vector or matrix")
- f_lapack.cgetrs_("N", &N, &NRHS, <float complex *>LU.data, &N,
- <l_int *>IPIV.data, <float complex *>b.data, &N,
- &info)
-
- assert info == 0, "Argument error in cgetrs"
-
- return b
-
-def zgetrs(np.ndarray[np.complex128_t, ndim=2] LU,
- np.ndarray[l_int, ndim=1] IPIV, B):
- cdef l_int N, NRHS, info
- cdef np.ndarray b
-
- assert_fortran_mat(LU)
-
- # again: workaround for 1x1-Fortran bug in NumPy < v2.0
- if (not isinstance(B, np.ndarray) or
- (B.ndim == 2 and (B.shape[0] > 1 or B.shape[1] > 1) and
- not B.flags["F_CONTIGUOUS"])):
- raise ValueError("In dgetrs: B must be a Fortran ordered NumPy array")
-
- b = B
- N = LU.shape[0]
- if b.ndim == 1:
+ if B.ndim == 1:
NRHS = 1
- elif b.ndim == 2:
- NRHS = b.shape[1]
- else:
- raise ValueError("In zgetrs: B must be a vector or matrix")
-
- f_lapack.zgetrs_("N", &N, &NRHS, <double complex *>LU.data, &N,
- <l_int *>IPIV.data, <double complex *>b.data, &N,
- &info)
-
- assert info == 0, "Argument error in zgetrs"
-
- return b
-
-# Wrappers for xGECON
+ elif B.ndim == 2:
+ NRHS = B.shape[1]
-def sgecon(np.ndarray[np.float32_t, ndim=2] LU,
- float normA, char *norm = b"1"):
+ if scalar is float:
+ lapack.sgetrs("N", &N, &NRHS, <float *>LU.data, &N,
+ <l_int *>IPIV.data, <float *>B.data, &N,
+ &info)
+ elif scalar is double:
+ lapack.dgetrs("N", &N, &NRHS, <double *>LU.data, &N,
+ <l_int *>IPIV.data, <double *>B.data, &N,
+ &info)
+ elif scalar is float_complex:
+ lapack.cgetrs("N", &N, &NRHS, <float complex *>LU.data, &N,
+ <l_int *>IPIV.data, <float complex *>B.data, &N,
+ &info)
+ elif scalar is double_complex:
+ lapack.zgetrs("N", &N, &NRHS, <double complex *>LU.data, &N,
+ <l_int *>IPIV.data, <double complex *>B.data, &N,
+ &info)
+
+ assert info == 0, "Argument error in getrs"
+
+ return B
+
+
+def gecon(np.ndarray[scalar, ndim=2] LU, double 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
+ cdef float srcond, snormA
+ cdef double drcond
- assert_fortran_mat(LU)
-
- N = LU.shape[0]
- work = np.empty(4*N, dtype = np.float32)
- iwork = np.empty(N, dtype = f_lapack.l_int_dtype)
-
- f_lapack.sgecon_(norm, &N, <float *>LU.data, &N, &normA,
- &rcond, <float *>work.data,
- <l_int *>iwork.data, &info)
-
- assert info == 0, "Argument error in sgecon"
-
- return rcond
-
-def dgecon(np.ndarray[np.float64_t, ndim=2] LU,
- 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
+ # Parameter checks
assert_fortran_mat(LU)
+ if norm[0] != b"1" and norm[0] != b"I":
+ raise ValueError("'norm' must be either '1' or 'I'")
+ if scalar in single_precision:
+ snormA = normA
- N = LU.shape[0]
- work = np.empty(4*N, dtype = np.float64)
- iwork = np.empty(N, dtype = f_lapack.l_int_dtype)
-
- f_lapack.dgecon_(norm, &N, <double *>LU.data, &N, &normA,
- &rcond, <double *>work.data,
- <l_int *>iwork.data, &info)
-
- assert info == 0, "Argument error in dgecon"
-
- return rcond
-
-def cgecon(np.ndarray[np.complex64_t, ndim=2] LU,
- 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
-
- assert_fortran_mat(LU)
+ # Allocate workspaces
N = LU.shape[0]
- work = np.empty(2*N, dtype = np.complex64)
- rwork = np.empty(2*N, dtype = np.float32)
- f_lapack.cgecon_(norm, &N, <float complex *>LU.data, &N, &normA,
- &rcond, <float complex *>work.data,
- <float *>rwork.data, &info)
-
- assert info == 0, "Argument error in cgecon"
-
- return rcond
+ cdef np.ndarray[l_int] iwork
+ if scalar in floating:
+ iwork = np.empty(N, dtype=int_dtype)
-def zgecon(np.ndarray[np.complex128_t, ndim=2] LU,
- 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
+ cdef np.ndarray[scalar] work
+ if scalar in floating:
+ work = np.empty(4 * N, dtype=LU.dtype)
+ else:
+ work = np.empty(2 * N, dtype=LU.dtype)
- assert_fortran_mat(LU)
+ cdef np.ndarray rwork
+ if scalar is float_complex:
+ rwork = np.empty(2 * N, dtype=np.float32)
+ elif scalar is double_complex:
+ rwork = np.empty(2 * N, dtype=np.float64)
- N = LU.shape[0]
- work = np.empty(2*N, dtype = np.complex128)
- rwork = np.empty(2*N, dtype = np.float64)
+ # The actual calculation
- f_lapack.zgecon_(norm, &N, <double complex *>LU.data, &N, &normA,
- &rcond, <double complex *>work.data,
- <double *>rwork.data, &info)
+ if scalar is float:
+ lapack.sgecon(norm, &N, <float *>LU.data, &N, &snormA,
+ &srcond, <float *>work.data,
+ <l_int *>iwork.data, &info)
+ elif scalar is double:
+ lapack.dgecon(norm, &N, <double *>LU.data, &N, &normA,
+ &drcond, <double *>work.data,
+ <l_int *>iwork.data, &info)
+ elif scalar is float_complex:
+ lapack.cgecon(norm, &N, <float complex *>LU.data, &N, &snormA,
+ &srcond, <float complex *>work.data,
+ <float *>rwork.data, &info)
+ elif scalar is double_complex:
+ lapack.zgecon(norm, &N, <double complex *>LU.data, &N, &normA,
+ &drcond, <double complex *>work.data,
+ <double *>rwork.data, &info)
- assert info == 0, "Argument error in zgecon"
+ assert info == 0, "Argument error in gecon"
- return rcond
+ if scalar in single_precision:
+ return srcond
+ else:
+ return drcond
-# Wrappers for xGGEV
# Helper function for xGGEV
def ggev_postprocess(dtype, alphar, alphai, vl_r=None, vr_r=None):
@@ -363,707 +264,363 @@ def ggev_postprocess(dtype, alphar, alphai, vl_r=None, vr_r=None):
return (alpha, vl, vr)
-def sggev(np.ndarray[np.float32_t, ndim=2] A,
- np.ndarray[np.float32_t, ndim=2] B,
- left=False, right=True):
+def ggev(np.ndarray[scalar, ndim=2] A, np.ndarray[scalar, 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
-
- assert_fortran_mat(A, B)
-
- N = A.shape[0]
- alphar = np.empty(N, dtype = np.float32)
- alphai = np.empty(N, dtype = np.float32)
- beta = np.empty(N, dtype = np.float32)
-
- if left:
- vl_r = np.empty((N,N), dtype = np.float32, order='F')
- vl_ptr = <float *>vl_r.data
- jobvl = "V"
- else:
- vl_r = None
- vl_ptr = NULL
- jobvl = "N"
-
- if right:
- vr_r = np.empty((N,N), dtype = np.float32, order='F')
- vr_ptr = <float *>vr_r.data
- jobvr = "V"
- else:
- vr_r = None
- vr_ptr = NULL
- jobvr = "N"
-
- # workspace query
- lwork = -1
-
- f_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)
- assert info == 0, "Argument error in sggev"
+ # Parameter checks
- lwork = <l_int>qwork
- work = np.empty(lwork, dtype = np.float32)
-
- # Now the real calculation
- f_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)
-
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in sggev")
-
- assert info == 0, "Argument error in sggev"
-
- alpha, vl, vr = ggev_postprocess(np.complex64, alphar, alphai, vl_r, vr_r)
+ assert_fortran_mat(A, B)
- return filter_args((True, True, left, right), (alpha, beta, vl, vr))
+ 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")
-def dggev(np.ndarray[np.float64_t, ndim=2] A,
- 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
+ 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")
- assert_fortran_mat(A, B)
+ # Allocate workspaces
N = A.shape[0]
- alphar = np.empty(N, dtype = np.float64)
- alphai = np.empty(N, dtype = np.float64)
- beta = np.empty(N, dtype = np.float64)
-
- if left:
- vl_r = np.empty((N,N), dtype = np.float64, order='F')
- vl_ptr = <double *>vl_r.data
- jobvl = "V"
- else:
- vl_r = None
- vl_ptr = NULL
- jobvl = "N"
- if right:
- vr_r = np.empty((N,N), dtype = np.float64, order='F')
- vr_ptr = <double *>vr_r.data
- jobvr = "V"
+ cdef np.ndarray[scalar] alphar, alphai
+ if scalar in cmplx:
+ alphar = np.empty(N, dtype=A.dtype)
+ alphai = None
else:
- vr_r = None
- vr_ptr = NULL
- jobvr = "N"
-
- # workspace query
- lwork = -1
-
- f_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)
-
- assert info == 0, "Argument error in dggev"
-
- lwork = <l_int>qwork
- work = np.empty(lwork, dtype = np.float64)
-
- # Now the real calculation
- f_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)
-
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in dggev")
-
- assert info == 0, "Argument error in dggev"
+ alphar = np.empty(N, dtype=A.dtype)
+ alphai = np.empty(N, dtype=A.dtype)
- alpha, vl, vr = ggev_postprocess(np.complex128, alphar, alphai, vl_r, vr_r)
-
- return filter_args((True, True, left, right), (alpha, beta, vl, vr))
+ 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)
-def cggev(np.ndarray[np.complex64_t, ndim=2] A,
- np.ndarray[np.complex64_t, ndim=2] B,
- left=False, right=True):
- cdef l_int N, info, lwork
+ cdef np.ndarray vl
+ cdef scalar *vl_ptr
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
-
- assert_fortran_mat(A, B)
-
- N = A.shape[0]
- alpha = np.empty(N, dtype = np.complex64)
- beta = np.empty(N, dtype = np.complex64)
-
if left:
- vl = np.empty((N,N), dtype = np.complex64, order='F')
- vl_ptr = <float complex *>vl.data
+ 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 = np.complex64, order='F')
- vr_ptr = <float complex *>vr.data
+ 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"
- rwork = np.empty(8 * N, dtype = np.float32)
-
- # workspace query
+ # Workspace query
+ # Xggev expects &qwork as a <scalar *> (even though it's an integer)
lwork = -1
- work = np.empty(1, dtype = np.complex64)
+ cdef scalar qwork
- f_lapack.cggev_(jobvl, jobvr, &N, <float complex *>A.data, &N,
- <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)
+ 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 cggev"
+ assert info == 0, "Argument error in ggev"
- lwork = <l_int>qwork.real
+ lwork = lwork_from_qwork(qwork)
+ cdef np.ndarray[scalar] work = np.empty(lwork, dtype=A.dtype)
- work = np.empty(lwork, dtype = np.complex64)
+ # The actual calculation
- # Now the real calculation
- f_lapack.cggev_(jobvl, jobvr, &N, <float complex *>A.data, &N,
- <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)
+ 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 cggev")
-
- assert info == 0, "Argument error in cggev"
-
- return filter_args((True, True, left, right), (alpha, beta, vl, vr))
-
-
-def zggev(np.ndarray[np.complex128_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A, B)
+ raise LinAlgError("QZ iteration failed to converge in sggev")
- N = A.shape[0]
- alpha = np.empty(N, dtype = np.complex128)
- beta = np.empty(N, dtype = np.complex128)
+ assert info == 0, "Argument error in ggev"
- if left:
- vl = np.empty((N,N), dtype = np.complex128, order='F')
- vl_ptr = <double complex *>vl.data
- jobvl = "V"
- else:
- vl_ptr = NULL
- jobvl = "N"
+ if scalar is float:
+ post_dtype = np.complex64
+ elif scalar is double:
+ post_dtype = np.complex128
- if right:
- vr = np.empty((N,N), dtype = np.complex128, order='F')
- vr_ptr = <double complex *>vr.data
- jobvr = "V"
- else:
- vr_ptr = NULL
- jobvr = "N"
+ cdef np.ndarray alpha
+ alpha = alphar
+ if scalar in floating:
+ alpha, vl, vr = ggev_postprocess(post_dtype, alphar, alphai, vl, vr)
- rwork = np.empty(8 * N, dtype = np.float64)
+ return filter_args((True, True, left, right), (alpha, beta, vl, vr))
- # workspace query
- lwork = -1
- work = np.empty(1, dtype = np.complex128)
- f_lapack.zggev_(jobvl, jobvr, &N, <double complex *>A.data, &N,
- <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)
+def gees(np.ndarray[scalar, ndim=2] A, calc_q=True, calc_ev=True):
+ cdef l_int N, lwork, sdim, info
- assert info == 0, "Argument error in zggev"
+ assert_fortran_mat(A)
- lwork = <l_int>qwork.real
- work = np.empty(lwork, dtype = np.complex128)
+ if A.ndim != 2:
+ raise ValueError("Expect matrix as input")
- # Now the real calculation
- f_lapack.zggev_(jobvl, jobvr, &N, <double complex *>A.data, &N,
- <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)
+ if A.shape[0] != A.shape[1]:
+ raise ValueError("Expect square matrix")
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in zggev")
+ # Allocate workspaces
- assert info == 0, "Argument error in zggev"
+ N = A.shape[0]
- return filter_args((True, True, left, right), (alpha, beta, vl, vr))
+ cdef np.ndarray[scalar] wr, wi
+ if scalar in cmplx:
+ wr = np.empty(N, dtype=A.dtype)
+ wi = None
+ else:
+ wr = np.empty(N, dtype=A.dtype)
+ wi = np.empty(N, dtype=A.dtype)
+ cdef np.ndarray rwork
+ if scalar is float_complex:
+ rwork = np.empty(N, dtype=np.float32)
+ elif scalar is double_complex:
+ rwork = np.empty(N, dtype=np.float64)
-# Wrapper for xGEES
-def sgees(np.ndarray[np.float32_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A)
-
- N = A.shape[0]
- wr = np.empty(N, dtype = np.float32)
- wi = np.empty(N, dtype = np.float32)
-
+ cdef scalar *vs_ptr
+ cdef np.ndarray[scalar, ndim=2] vs
if calc_q:
- vs = np.empty((N,N), dtype = np.float32, order='F')
- vs_ptr = <float *>vs.data
+ vs = np.empty((N,N), dtype=A.dtype, order='F')
+ vs_ptr = <scalar *>vs.data
jobvs = "V"
else:
+ vs = None
vs_ptr = NULL
jobvs = "N"
- # workspace query
+ # Workspace query
+ # Xgees expects &qwork as a <scalar *> (even though it's an integer)
lwork = -1
- f_lapack.sgees_(jobvs, "N", NULL, &N, <float *>A.data, &N,
- &sdim, <float *>wr.data, <float *>wi.data, vs_ptr, &N,
- &qwork, &lwork, NULL, &info)
+ cdef scalar qwork
+
+ if scalar is float:
+ lapack.sgees(jobvs, "N", NULL, &N, <float *>A.data, &N,
+ &sdim, <float *>wr.data, <float *>wi.data, vs_ptr, &N,
+ &qwork, &lwork, NULL, &info)
+ elif scalar is double:
+ lapack.dgees(jobvs, "N", NULL, &N, <double *>A.data, &N,
+ &sdim, <double *>wr.data, <double *>wi.data, vs_ptr, &N,
+ &qwork, &lwork, NULL, &info)
+ elif scalar is float_complex:
+ lapack.cgees(jobvs, "N", NULL, &N, <float complex *>A.data, &N,
+ &sdim, <float complex *>wr.data, vs_ptr, &N,
+ &qwork, &lwork, <float *>rwork.data, NULL, &info)
+ elif scalar is double_complex:
+ lapack.zgees(jobvs, "N", NULL, &N, <double complex *>A.data, &N,
+ &sdim, <double complex *>wr.data, vs_ptr, &N,
+ &qwork, &lwork, <double *>rwork.data, NULL, &info)
assert info == 0, "Argument error in sgees"
- lwork = <int>qwork
- work = np.empty(lwork, dtype = np.float32)
-
- # Now the real calculation
- f_lapack.sgees_(jobvs, "N", NULL, &N, <float *>A.data, &N,
- &sdim, <float *>wr.data, <float *>wi.data, vs_ptr, &N,
- <float *>work.data, &lwork, NULL, &info)
+ lwork = lwork_from_qwork(qwork)
+ cdef np.ndarray[scalar] work = np.empty(lwork, dtype=A.dtype)
+
+ # The actual calculation
+
+ if scalar is float:
+ lapack.sgees(jobvs, "N", NULL, &N, <float *>A.data, &N,
+ &sdim, <float *>wr.data, <float *>wi.data, vs_ptr, &N,
+ <float *>work.data, &lwork, NULL, &info)
+ elif scalar is double:
+ lapack.dgees(jobvs, "N", NULL, &N, <double *>A.data, &N,
+ &sdim, <double *>wr.data, <double *>wi.data, vs_ptr, &N,
+ <double *>work.data, &lwork, NULL, &info)
+ elif scalar is float_complex:
+ lapack.cgees(jobvs, "N", NULL, &N, <float complex *>A.data, &N,
+ &sdim, <float complex *>wr.data, vs_ptr, &N,
+ <float complex *>work.data, &lwork,
+ <float *>rwork.data, NULL, &info)
+ elif scalar is double_complex:
+ lapack.zgees(jobvs, "N", NULL, &N, <double complex *>A.data, &N,
+ &sdim, <double complex *>wr.data, vs_ptr, &N,
+ <double complex *>work.data, &lwork,
+ <double *>rwork.data, NULL, &info)
if info > 0:
- raise LinAlgError("QR iteration failed to converge in sgees")
+ raise LinAlgError("QR iteration failed to converge in gees")
- assert info == 0, "Argument error in sgees"
+ assert info == 0, "Argument error in gees"
- if wi.nonzero()[0].size:
- w = wr + 1j * wi
- else:
- w = wr
+ # Real inputs possibly produce complex output
+ cdef np.ndarray w = maybe_complex[scalar](0, wr, wi)
return filter_args((True, calc_q, calc_ev), (A, vs, w))
-def dgees(np.ndarray[np.float64_t, ndim=2] A,
- 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
+def trsen(np.ndarray[l_logical] select,
+ np.ndarray[scalar, ndim=2] T,
+ np.ndarray[scalar, ndim=2] Q,
+ calc_ev=True):
+ cdef l_int N, M, lwork, liwork, qiwork, info
- assert_fortran_mat(A)
+ assert_fortran_mat(T, Q)
- N = A.shape[0]
- wr = np.empty(N, dtype = np.float64)
- wi = np.empty(N, dtype = np.float64)
+ # Allocate workspaces
- if calc_q:
- vs = np.empty((N,N), dtype = np.float64, order='F')
- vs_ptr = <double *>vs.data
- jobvs = "V"
- else:
- vs_ptr = NULL
- jobvs = "N"
+ N = T.shape[0]
- # workspace query
- lwork = -1
- f_lapack.dgees_(jobvs, "N", NULL, &N, <double *>A.data, &N,
- &sdim, <double *>wr.data, <double *>wi.data, vs_ptr, &N,
- &qwork, &lwork, NULL, &info)
+ cdef np.ndarray[scalar] wr, wi
+ if scalar in cmplx:
+ wr = np.empty(N, dtype=T.dtype)
+ wi = None
+ else:
+ wr = np.empty(N, dtype=T.dtype)
+ wi = np.empty(N, dtype=T.dtype)
- assert info == 0, "Argument error in dgees"
+ cdef char *compq
+ cdef scalar *q_ptr
+ if Q is not None:
+ compq = "V"
+ q_ptr = <scalar *>Q.data
+ else:
+ compq = "N"
+ q_ptr = NULL
- lwork = <int>qwork
- work = np.empty(lwork, dtype = np.float64)
+ # Workspace query
+ # Xtrsen expects &qwork as a <scalar *> (even though it's an integer)
+ cdef scalar qwork
+ lwork = liwork = -1
- # Now the real calculation
- f_lapack.dgees_(jobvs, "N", NULL, &N, <double *>A.data, &N,
- &sdim, <double *>wr.data, <double *>wi.data, vs_ptr, &N,
- <double *>work.data, &lwork, NULL, &info)
+ if scalar is float:
+ 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)
+ elif scalar is double:
+ lapack.dtrsen("N", compq, <l_logical *>select.data,
+ &N, <double *>T.data, &N, q_ptr, &N,
+ <double *>wr.data, <double *>wi.data, &M, NULL, NULL,
+ &qwork, &lwork, &qiwork, &liwork, &info)
+ elif scalar is float_complex:
+ lapack.ctrsen("N", compq, <l_logical *>select.data,
+ &N, <float complex *>T.data, &N, q_ptr, &N,
+ <float complex *>wr.data, &M, NULL, NULL,
+ &qwork, &lwork, &info)
+ elif scalar is double_complex:
+ lapack.ztrsen("N", compq, <l_logical *>select.data,
+ &N, <double complex *>T.data, &N, q_ptr, &N,
+ <double complex *>wr.data, &M, NULL, NULL,
+ &qwork, &lwork, &info)
+
+ assert info == 0, "Argument error in trsen"
+
+ lwork = lwork_from_qwork(qwork)
+ cdef np.ndarray[scalar] work = np.empty(lwork, dtype=T.dtype)
+
+ cdef np.ndarray[l_int] iwork = None
+ if scalar in floating:
+ liwork = qiwork
+ iwork = np.empty(liwork, dtype=int_dtype)
+
+ # Tha actual calculation
+
+ if scalar is float:
+ 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,
+ <l_int *>iwork.data, &liwork, &info)
+ elif scalar is double:
+ lapack.dtrsen("N", compq, <l_logical *>select.data,
+ &N, <double *>T.data, &N, q_ptr, &N,
+ <double *>wr.data, <double *>wi.data, &M, NULL, NULL,
+ <double *>work.data, &lwork,
+ <l_int *>iwork.data, &liwork, &info)
+ elif scalar is float_complex:
+ lapack.ctrsen("N", compq, <l_logical *>select.data,
+ &N, <float complex *>T.data, &N, q_ptr, &N,
+ <float complex *>wr.data, &M, NULL, NULL,
+ <float complex *>work.data, &lwork, &info)
+ elif scalar is double_complex:
+ lapack.ztrsen("N", compq, <l_logical *>select.data,
+ &N, <double complex *>T.data, &N, q_ptr, &N,
+ <double complex *>wr.data, &M, NULL, NULL,
+ <double complex *>work.data, &lwork, &info)
if info > 0:
- raise LinAlgError("QR iteration failed to converge in dgees")
+ raise LinAlgError("Reordering failed; problem is very ill-conditioned")
- assert info == 0, "Argument error in dgees"
+ assert info == 0, "Argument error in trsen"
- if wi.nonzero()[0].size:
- w = wr + 1j * wi
- else:
- w = wr
+ # Real inputs possibly produce complex output
+ cdef np.ndarray w = maybe_complex[scalar](0, wr, wi)
- return filter_args((True, calc_q, calc_ev), (A, vs, w))
+ return filter_args((True, Q is not None, calc_ev), (T, Q, w))
-def cgees(np.ndarray[np.complex64_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A)
-
- N = A.shape[0]
- w = np.empty(N, dtype = np.complex64)
- rwork = np.empty(N, dtype = np.float32)
-
- if calc_q:
- vs = np.empty((N,N), dtype = np.complex64, order='F')
- vs_ptr = <float complex *>vs.data
- jobvs = "V"
- else:
- vs_ptr = NULL
- jobvs = "N"
-
- # workspace query
- lwork = -1
- f_lapack.cgees_(jobvs, "N", NULL, &N, <float complex *>A.data, &N,
- &sdim, <float complex *>w.data, vs_ptr, &N,
- &qwork, &lwork, <float *>rwork.data, NULL, &info)
-
- assert info == 0, "Argument error in cgees"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex64)
-
- # Now the real calculation
- f_lapack.cgees_(jobvs, "N", NULL, &N, <float complex *>A.data, &N,
- &sdim, <float complex *>w.data, vs_ptr, &N,
- <float complex *>work.data, &lwork,
- <float *>rwork.data, NULL, &info)
-
- if info > 0:
- raise LinAlgError("QR iteration failed to converge in cgees")
-
- assert info == 0, "Argument error in cgees"
-
- return filter_args((True, calc_q, calc_ev), (A, vs, w))
-
-
-def zgees(np.ndarray[np.complex128_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A)
-
- N = A.shape[0]
- w = np.empty(N, dtype = np.complex128)
- rwork = np.empty(N, dtype = np.float64)
-
- if calc_q:
- vs = np.empty((N,N), dtype = np.complex128, order='F')
- vs_ptr = <double complex *>vs.data
- jobvs = "V"
- else:
- vs_ptr = NULL
- jobvs = "N"
-
- # workspace query
- lwork = -1
- f_lapack.zgees_(jobvs, "N", NULL, &N, <double complex *>A.data, &N,
- &sdim, <double complex *>w.data, vs_ptr, &N,
- &qwork, &lwork, <double *>rwork.data, NULL, &info)
-
- assert info == 0, "Argument error in zgees"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex128)
-
- # Now the real calculation
- f_lapack.zgees_(jobvs, "N", NULL, &N, <double complex *>A.data, &N,
- &sdim, <double complex *>w.data, vs_ptr, &N,
- <double complex *>work.data, &lwork,
- <double *>rwork.data, NULL, &info)
-
- if info > 0:
- raise LinAlgError("QR iteration failed to converge in zgees")
-
- assert info == 0, "Argument error in zgees"
-
- return filter_args((True, calc_q, calc_ev), (A, vs, w))
-
-# Wrapper for xTRSEN
-def strsen(np.ndarray[l_logical] select,
- 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))
-
-
-def dtrsen(np.ndarray[l_logical] select,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- wr = np.empty(N, dtype = np.float64)
- wi = np.empty(N, dtype = np.float64)
-
- if Q is not None:
- compq = "V"
- q_ptr = <double *>Q.data
- else:
- compq = "N"
- q_ptr = NULL
-
- # workspace query
- lwork = liwork = -1
- f_lapack.dtrsen_("N", compq, <l_logical *>select.data,
- &N, <double *>T.data, &N, q_ptr, &N,
- <double *>wr.data, <double *>wi.data, &M, NULL, NULL,
- &qwork, &lwork, &qiwork, &liwork, &info)
-
- assert info == 0, "Argument error in dtrsen"
-
- lwork = <int>qwork
- work = np.empty(lwork, dtype = np.float64)
- liwork = qiwork
- iwork = np.empty(liwork, dtype = f_lapack.l_int_dtype)
-
- # Now the real calculation
- f_lapack.dtrsen_("N", compq, <l_logical *>select.data,
- &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)
-
- if info > 0:
- raise LinAlgError("Reordering failed; problem is very ill-conditioned")
-
- assert info == 0, "Argument error in dtrsen"
-
- 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))
-
-
-def ctrsen(np.ndarray[l_logical] select,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- w = np.empty(N, dtype = np.complex64)
-
- if Q is not None:
- compq = "V"
- q_ptr = <float complex *>Q.data
- else:
- compq = "N"
- q_ptr = NULL
-
- # workspace query
- lwork = -1
- f_lapack.ctrsen_("N", compq, <l_logical *>select.data,
- &N, <float complex *>T.data, &N, q_ptr, &N,
- <float complex *>w.data, &M, NULL, NULL,
- &qwork, &lwork, &info)
-
- assert info == 0, "Argument error in ctrsen"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex64)
-
- # Now the real calculation
- f_lapack.ctrsen_("N", compq, <l_logical *>select.data,
- &N, <float complex *>T.data, &N, q_ptr, &N,
- <float complex *>w.data, &M, NULL, NULL,
- <float complex *>work.data, &lwork, &info)
-
- if info > 0:
- raise LinAlgError("Reordering failed; problem is very ill-conditioned")
-
- assert info == 0, "Argument error in ctrsen"
-
- return filter_args((True, Q is not None, calc_ev), (T, Q, w))
-
-
-def ztrsen(np.ndarray[l_logical] select,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- w = np.empty(N, dtype = np.complex128)
-
- if Q is not None:
- compq = "V"
- q_ptr = <double complex *>Q.data
- else:
- compq = "N"
- q_ptr = NULL
-
- # workspace query
- lwork = -1
- f_lapack.ztrsen_("N", compq, <l_logical *>select.data,
- &N, <double complex *>T.data, &N, q_ptr, &N,
- <double complex *>w.data, &M, NULL, NULL,
- &qwork, &lwork, &info)
-
- assert info == 0, "Argument error in ztrsen"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex128)
-
- # Now the real calculation
- f_lapack.ztrsen_("N", compq, <l_logical *>select.data,
- &N, <double complex *>T.data, &N, q_ptr, &N,
- <double complex *>w.data, &M, NULL, NULL,
- <double complex *>work.data, &lwork, &info)
-
- if info > 0:
- raise LinAlgError("Reordering failed; problem is very ill-conditioned")
-
- assert info == 0, "Argument error in ztrsen"
-
- return filter_args((True, Q is not None, calc_ev), (T, Q, w))
-
-
-# Helper function for xTREVC and xTGEVC
-def txevc_postprocess(dtype, T, vreal, np.ndarray[l_logical] select):
- cdef int N, M, i, m, indx
+# Helper function for xTREVC and xTGEVC
+def txevc_postprocess(dtype, T, vreal, np.ndarray[l_logical] select):
+ cdef int N, M, i, m, indx
N = T.shape[0]
if select is None:
- select = np.ones(N, dtype = f_lapack.l_logical_dtype)
+ select = np.ones(N, dtype = logical_dtype)
selindx = select.nonzero()[0]
M = selindx.size
@@ -1085,1000 +642,48 @@ def txevc_postprocess(dtype, T, vreal, np.ndarray[l_logical] select):
elif k > 0 and T[k,k-1]:
# we have the situation of a 2x2 block, and
# the eigenvalue with the negative imaginary part desired
- v[:, m] = vreal[:, indx] - 1j * vreal[:, indx + 1]
-
- indx += 2
- else:
- # real eigenvalue
- v[:, m] = vreal[:, indx]
-
- indx += 1
- return v
-
-
-# Wrappers for xTREVC
-def strevc(np.ndarray[np.float32_t, ndim=2] T,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- work = np.empty(4*N, dtype = np.float32)
-
- if left and right:
- side = "B"
- elif left:
- side = "L"
- elif right:
- side = "R"
- else:
- return
-
- if select is not None:
- howmny = "S"
- MM = select.nonzero()[0].size
- # Correct for possible additional storage if a single complex
- # eigenvalue is selected.
- # For that: Figure out the positions of the 2x2 blocks.
- cmplxindx = np.diagonal(T, -1).nonzero()[0]
- for i in cmplxindx:
- if bool(select[i]) != bool(select[i+1]):
- MM += 1
-
- # Select is overwritten in strevc.
- select_cpy = np.array(select, dtype = f_lapack.l_logical_dtype,
- order = 'F')
- select_ptr = <l_logical *>select_cpy.data
- else:
- MM = N
- select_ptr = NULL
- if Q is not None:
- howmny = "B"
- else:
- howmny = "A"
-
- if left:
- if Q is not None and select is None:
- vl_r = np.asfortranarray(Q.copy())
- else:
- vl_r = np.empty((N, MM), dtype = np.float32, order='F')
- vl_r_ptr = <float *>vl_r.data
- else:
- vl_r_ptr = NULL
-
- if right:
- if Q is not None and select is None:
- vr_r = np.asfortranarray(Q.copy())
- else:
- vr_r = np.empty((N, MM), dtype = np.float32, order='F')
- vr_r_ptr = <float *>vr_r.data
- else:
- vr_r_ptr = NULL
-
- f_lapack.strevc_(side, howmny, select_ptr,
- &N, <float *>T.data, &N,
- vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
- <float *>work.data, &info)
-
- assert info == 0, "Argument error in strevc"
- assert MM == M, "Unexpected number of eigenvectors returned in strevc"
-
- if select is not None and Q is not None:
- if left:
- vl_r = np.asfortranarray(np.dot(Q, vl_r))
- if right:
- vr_r = np.asfortranarray(np.dot(Q, vr_r))
-
- # If there are complex eigenvalues, we need to postprocess the
- # eigenvectors.
- if np.diagonal(T, -1).nonzero()[0].size:
- if left:
- vl = txevc_postprocess(np.complex64, T, vl_r, select)
- if right:
- vr = txevc_postprocess(np.complex64, T, vr_r, select)
- else:
- if left:
- vl = vl_r
- if right:
- vr = vr_r
-
- if left and right:
- return (vl, vr)
- elif left:
- return vl
- else:
- return vr
-
-
-def dtrevc(np.ndarray[np.float64_t, ndim=2] T,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- work = np.empty(4*N, dtype = np.float64)
-
- if left and right:
- side = "B"
- elif left:
- side = "L"
- elif right:
- side = "R"
- else:
- return
-
- if select is not None:
- howmny = "S"
- MM = select.nonzero()[0].size
- # Correct for possible additional storage if a single complex
- # eigenvalue is selected.
- # For that: Figure out the positions of the 2x2 blocks.
- cmplxindx = np.diagonal(T, -1).nonzero()[0]
- for i in cmplxindx:
- if bool(select[i]) != bool(select[i+1]):
- MM += 1
-
- # Select is overwritten in dtrevc.
- select_cpy = np.array(select, dtype = f_lapack.l_logical_dtype,
- order = 'F')
- select_ptr = <l_logical *>select_cpy.data
- else:
- MM = N
- select_ptr = NULL
- if Q is not None:
- howmny = "B"
- else:
- howmny = "A"
-
- if left:
- if Q is not None and select is None:
- vl_r = np.asfortranarray(Q.copy())
- else:
- vl_r = np.empty((N, MM), dtype = np.float64, order='F')
- vl_r_ptr = <double *>vl_r.data
- else:
- vl_r_ptr = NULL
-
- if right:
- if Q is not None and select is None:
- vr_r = np.asfortranarray(Q.copy())
- else:
- vr_r = np.empty((N, MM), dtype = np.float64, order='F')
- vr_r_ptr = <double *>vr_r.data
- else:
- vr_r_ptr = NULL
-
- f_lapack.dtrevc_(side, howmny, select_ptr,
- &N, <double *>T.data, &N,
- vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
- <double *>work.data, &info)
-
- assert info == 0, "Argument error in dtrevc"
- assert MM == M, "Unexpected number of eigenvectors returned in dtrevc"
-
- if select is not None and Q is not None:
- if left:
- vl_r = np.asfortranarray(np.dot(Q, vl_r))
- if right:
- vr_r = np.asfortranarray(np.dot(Q, vr_r))
-
- # If there are complex eigenvalues, we need to postprocess the eigenvectors
- if np.diagonal(T, -1).nonzero()[0].size:
- if left:
- vl = txevc_postprocess(np.complex128, T, vl_r, select)
- if right:
- vr = txevc_postprocess(np.complex128, T, vr_r, select)
- else:
- if left:
- vl = vl_r
- if right:
- vr = vr_r
-
- if left and right:
- return (vl, vr)
- elif left:
- return vl
- else:
- return vr
-
-
-def ctrevc(np.ndarray[np.complex64_t, ndim=2] T,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- work = np.empty(2*N, dtype = np.complex64)
- rwork = np.empty(N, dtype = np.float32)
-
- if left and right:
- side = "B"
- elif left:
- side = "L"
- elif right:
- side = "R"
- else:
- return
-
- if select is not None:
- howmny = "S"
- MM = select.nonzero()[0].size
- select_ptr = <l_logical *>select.data
- else:
- MM = N
- select_ptr = NULL
- if Q is not None:
- howmny = "B"
- else:
- howmny = "A"
-
- if left:
- if Q is not None and select is None:
- vl = np.asfortranarray(Q.copy())
- else:
- vl = np.empty((N, MM), dtype = np.complex64, order='F')
- vl_ptr = <float complex *>vl.data
- else:
- vl_ptr = NULL
-
- if right:
- if Q is not None and select is None:
- vr = np.asfortranarray(Q.copy())
- else:
- vr = np.empty((N, MM), dtype = np.complex64, order='F')
- vr_ptr = <float complex *>vr.data
- else:
- vr_ptr = NULL
-
- f_lapack.ctrevc_(side, howmny, select_ptr,
- &N, <float complex *>T.data, &N,
- vl_ptr, &N, vr_ptr, &N, &MM, &M,
- <float complex *>work.data, <float *>rwork.data, &info)
-
- assert info == 0, "Argument error in ctrevc"
- assert MM == M, "Unexpected number of eigenvectors returned in ctrevc"
-
- if select is not None and Q is not None:
- if left:
- vl = np.asfortranarray(np.dot(Q, vl))
- if right:
- vr = np.asfortranarray(np.dot(Q, vr))
-
- if left and right:
- return (vl, vr)
- elif left:
- return vl
- else:
- return vr
-
-
-def ztrevc(np.ndarray[np.complex128_t, ndim=2] T,
- 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
-
- assert_fortran_mat(T, Q)
-
- N = T.shape[0]
- work = np.empty(2*N, dtype = np.complex128)
- rwork = np.empty(N, dtype = np.float64)
-
- if left and right:
- side = "B"
- elif left:
- side = "L"
- elif right:
- side = "R"
- else:
- return
-
- if select is not None:
- howmny = "S"
- MM = select.nonzero()[0].size
- select_ptr = <l_logical *>select.data
- else:
- MM = N
- select_ptr = NULL
- if Q is not None:
- howmny = "B"
- else:
- howmny = "A"
-
- if left:
- if Q is not None and select is None:
- vl = np.asfortranarray(Q.copy())
- else:
- vl = np.empty((N, MM), dtype = np.complex128, order='F')
- vl_ptr = <double complex *>vl.data
- else:
- vl_ptr = NULL
-
- if right:
- if Q is not None and select is None:
- vr = np.asfortranarray(Q.copy())
- else:
- vr = np.empty((N, MM), dtype = np.complex128, order='F')
- vr_ptr = <double complex *>vr.data
- else:
- vr_ptr = NULL
-
- f_lapack.ztrevc_(side, howmny, select_ptr,
- &N, <double complex *>T.data, &N,
- vl_ptr, &N, vr_ptr, &N, &MM, &M,
- <double complex *>work.data, <double *>rwork.data, &info)
-
- assert info == 0, "Argument error in ztrevc"
- assert MM == M, "Unexpected number of eigenvectors returned in ztrevc"
-
- if select is not None and Q is not None:
- if left:
- vl = np.asfortranarray(np.dot(Q, vl))
- if right:
- vr = np.asfortranarray(np.dot(Q, vr))
-
- if left and right:
- return (vl, vr)
- elif left:
- return vl
- else:
- return vr
-
-
-# wrappers for xGGES
-def sgges(np.ndarray[np.float32_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A, B)
-
- N = A.shape[0]
- alphar = np.empty(N, dtype = np.float32)
- alphai = np.empty(N, dtype = np.float32)
- beta = np.empty(N, dtype = np.float32)
-
- if calc_q:
- vsl = np.empty((N,N), dtype = np.float32, order='F')
- vsl_ptr = <float *>vsl.data
- jobvsl = "V"
- else:
- vsl = None
- vsl_ptr = NULL
- jobvsl = "N"
-
- if calc_z:
- vsr = np.empty((N,N), dtype = np.float32, order='F')
- vsr_ptr = <float *>vsr.data
- jobvsr = "V"
- else:
- vsr = None
- vsr_ptr = NULL
- jobvsr = "N"
-
- # workspace query
- lwork = -1
- f_lapack.sgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- assert info == 0, "Argument error in zgees"
-
- lwork = <int>qwork
- work = np.empty(lwork, dtype = np.float32)
-
- # Now the real calculation
- f_lapack.sgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in sgges")
-
- assert info == 0, "Argument error in zgees"
-
- if alphai.nonzero()[0].size:
- alpha = alphar + 1j * alphai
- else:
- alpha = alphar
-
- return filter_args((True, True, calc_q, calc_z, calc_ev, calc_ev),
- (A, B, vsl, vsr, alpha, beta))
-
-
-def dgges(np.ndarray[np.float64_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A, B)
-
- N = A.shape[0]
- alphar = np.empty(N, dtype = np.float64)
- alphai = np.empty(N, dtype = np.float64)
- beta = np.empty(N, dtype = np.float64)
-
- if calc_q:
- vsl = np.empty((N,N), dtype = np.float64, order='F')
- vsl_ptr = <double *>vsl.data
- jobvsl = "V"
- else:
- vsl = None
- vsl_ptr = NULL
- jobvsl = "N"
-
- if calc_z:
- vsr = np.empty((N,N), dtype = np.float64, order='F')
- vsr_ptr = <double *>vsr.data
- jobvsr = "V"
- else:
- vsr = None
- vsr_ptr = NULL
- jobvsr = "N"
-
- # workspace query
- lwork = -1
- f_lapack.dgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- assert info == 0, "Argument error in zgees"
-
- lwork = <int>qwork
- work = np.empty(lwork, dtype = np.float64)
-
- # Now the real calculation
- f_lapack.dgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in dgges")
-
- assert info == 0, "Argument error in zgees"
-
- if alphai.nonzero()[0].size:
- alpha = alphar + 1j * alphai
- else:
- alpha = alphar
-
- return filter_args((True, True, calc_q, calc_z, calc_ev, calc_ev),
- (A, B, vsl, vsr, alpha, beta))
-
-
-def cgges(np.ndarray[np.complex64_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A, B)
-
- N = A.shape[0]
- alpha = np.empty(N, dtype = np.complex64)
- beta = np.empty(N, dtype = np.complex64)
- rwork = np.empty(8*N, dtype = np.float32)
-
- if calc_q:
- vsl = np.empty((N,N), dtype = np.complex64, order='F')
- vsl_ptr = <float complex *>vsl.data
- jobvsl = "V"
- else:
- vsl = None
- vsl_ptr = NULL
- jobvsl = "N"
-
- if calc_z:
- vsr = np.empty((N,N), dtype = np.complex64, order='F')
- vsr_ptr = <float complex *>vsr.data
- jobvsr = "V"
- else:
- vsr = None
- vsr_ptr = NULL
- jobvsr = "N"
-
- # workspace query
- lwork = -1
- f_lapack.cgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- assert info == 0, "Argument error in zgees"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex64)
-
- # Now the real calculation
- f_lapack.cgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in cgges")
-
- assert info == 0, "Argument error in zgees"
-
- return filter_args((True, True, calc_q, calc_z, calc_ev, calc_ev),
- (A, B, vsl, vsr, alpha, beta))
-
-
-def zgges(np.ndarray[np.complex128_t, ndim=2] A,
- 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
-
- assert_fortran_mat(A, B)
-
- N = A.shape[0]
- alpha = np.empty(N, dtype = np.complex128)
- beta = np.empty(N, dtype = np.complex128)
- rwork = np.empty(8*N, dtype = np.float64)
-
- if calc_q:
- vsl = np.empty((N,N), dtype = np.complex128, order='F')
- vsl_ptr = <double complex *>vsl.data
- jobvsl = "V"
- else:
- vsl = None
- vsl_ptr = NULL
- jobvsl = "N"
-
- if calc_z:
- vsr = np.empty((N,N), dtype = np.complex128, order='F')
- vsr_ptr = <double complex *>vsr.data
- jobvsr = "V"
- else:
- vsr = None
- vsr_ptr = NULL
- jobvsr = "N"
-
- # workspace query
- lwork = -1
- f_lapack.zgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- assert info == 0, "Argument error in zgees"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex128)
-
- # Now the real calculation
- f_lapack.zgges_(jobvsl, jobvsr, "N", NULL,
- &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)
-
- if info > 0:
- raise LinAlgError("QZ iteration failed to converge in zgges")
-
- assert info == 0, "Argument error in zgees"
-
- return filter_args((True, True, calc_q, calc_z, calc_ev, calc_ev),
- (A, B, vsl, vsr, alpha, beta))
-
-
-# wrappers for xTGSEN
-def stgsen(np.ndarray[l_logical] select,
- 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))
-
-
-def dtgsen(np.ndarray[l_logical] select,
- 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
-
- assert_fortran_mat(S, T, Q, Z)
-
- N = S.shape[0]
- alphar = np.empty(N, dtype = np.float64)
- alphai = np.empty(N, dtype = np.float64)
- beta = np.empty(N, dtype = np.float64)
- ijob = 0
-
- if Q is not None:
- wantq = 1
- q_ptr = <double *>Q.data
- else:
- wantq = 0
- q_ptr = NULL
-
- if Z is not None:
- wantz = 1
- z_ptr = <double *>Z.data
- else:
- wantz = 0
- z_ptr = NULL
-
- # workspace query
- lwork = -1
- liwork = -1
- f_lapack.dtgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data,
- &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)
-
- assert info == 0, "Argument error in dtgsen"
-
- lwork = <int>qwork
- work = np.empty(lwork, dtype = np.float64)
- liwork = qiwork
- iwork = np.empty(liwork, dtype = int_dtype)
-
- # Now the real calculation
- f_lapack.dtgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data,
- &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)
-
- if info > 0:
- raise LinAlgError("Reordering failed; problem is very ill-conditioned")
-
- assert info == 0, "Argument error in dtgsen"
-
- 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))
-
-
-def ctgsen(np.ndarray[l_logical] select,
- 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
-
- assert_fortran_mat(S, T, Q, Z)
-
- N = S.shape[0]
- alpha = np.empty(N, dtype = np.complex64)
- beta = np.empty(N, dtype = np.complex64)
- ijob = 0
-
- if Q is not None:
- wantq = 1
- q_ptr = <float complex *>Q.data
- else:
- wantq = 0
- q_ptr = NULL
-
- if Z is not None:
- wantz = 1
- z_ptr = <float complex *>Z.data
- else:
- wantz = 0
- z_ptr = NULL
-
- # workspace query
- lwork = -1
- liwork = -1
- f_lapack.ctgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data,
- &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)
-
- assert info == 0, "Argument error in ctgsen"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex64)
- liwork = qiwork
- iwork = np.empty(liwork, dtype = int_dtype)
-
- # Now the real calculation
- f_lapack.ctgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data,
- &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)
-
- if info > 0:
- raise LinAlgError("Reordering failed; problem is very ill-conditioned")
-
- assert info == 0, "Argument error in ctgsen"
-
- return filter_args((True, True, Q is not None, Z is not None,
- calc_ev, calc_ev),
- (S, T, Q, Z, alpha, beta))
-
-
-def ztgsen(np.ndarray[l_logical] select,
- 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
-
- assert_fortran_mat(S, T, Q, Z)
-
- N = S.shape[0]
- alpha = np.empty(N, dtype = np.complex128)
- beta = np.empty(N, dtype = np.complex128)
- ijob = 0
-
- if Q is not None:
- wantq = 1
- q_ptr = <double complex *>Q.data
- else:
- wantq = 0
- q_ptr = NULL
-
- if Z is not None:
- wantz = 1
- z_ptr = <double complex *>Z.data
- else:
- wantz = 0
- z_ptr = NULL
-
- # workspace query
- lwork = -1
- liwork = -1
- f_lapack.ztgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data,
- &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)
-
- assert info == 0, "Argument error in ztgsen"
-
- lwork = <int>qwork.real
- work = np.empty(lwork, dtype = np.complex128)
- liwork = qiwork
- iwork = np.empty(liwork, dtype = int_dtype)
-
- # Now the real calculation
- f_lapack.ztgsen_(&ijob, &wantq, &wantz, <l_logical *>select.data,
- &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)
-
- if info > 0:
- raise LinAlgError("Reordering failed; problem is very ill-conditioned")
+ v[:, m] = vreal[:, indx] - 1j * vreal[:, indx + 1]
- assert info == 0, "Argument error in ztgsen"
+ indx += 2
+ else:
+ # real eigenvalue
+ v[:, m] = vreal[:, indx]
- return filter_args((True, True, Q is not None, Z is not None,
- calc_ev, calc_ev),
- (S, T, Q, Z, alpha, beta))
+ indx += 1
+ return v
-# xTGEVC
-def stgevc(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,
- np.ndarray[l_logical] select=None,
- left=False, right=True):
+def trevc(np.ndarray[scalar, ndim=2] T,
+ np.ndarray[scalar, ndim=2] Q,
+ np.ndarray[l_logical] select,
+ 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
- assert_fortran_mat(S, T, Q, Z)
+ # Parameter checks
- N = S.shape[0]
- work = np.empty(6*N, dtype = np.float32)
+ if (T.shape[0] != T.shape[1] or Q.shape[0] != Q.shape[1]
+ or T.shape[0] != Q.shape[0]):
+ raise ValueError("Invalid Schur decomposition as input")
+
+ assert_fortran_mat(T, Q)
+
+ # Workspace allocation
+
+ N = T.shape[0]
+
+ cdef np.ndarray[scalar] work
+ if scalar in floating:
+ work = np.empty(4 * N, dtype=T.dtype)
+ else:
+ work = np.empty(2 * N, dtype=T.dtype)
+
+ cdef np.ndarray rwork = None
+ if scalar is float_complex:
+ rwork = np.empty(N, dtype=np.float32)
+ elif scalar is double_complex:
+ rwork = np.empty(N, dtype=np.float64)
if left and right:
side = "B"
@@ -2089,78 +694,102 @@ def stgevc(np.ndarray[np.float32_t, ndim=2] S,
else:
return
- backtr = False
-
+ cdef np.ndarray[l_logical] select_cpy
+ cdef l_logical *select_ptr
if select is not None:
howmny = "S"
MM = select.nonzero()[0].size
# Correct for possible additional storage if a single complex
# eigenvalue is selected.
# For that: Figure out the positions of the 2x2 blocks.
- cmplxindx = np.diagonal(S, -1).nonzero()[0]
+ cmplxindx = np.diagonal(T, -1).nonzero()[0]
for i in cmplxindx:
if bool(select[i]) != bool(select[i+1]):
MM += 1
- # select is overwritten in stgevc
- select_cpy = np.array(select, dtype = f_lapack.l_logical_dtype,
+ # Select is overwritten in strevc.
+ select_cpy = np.array(select, dtype = logical_dtype,
order = 'F')
select_ptr = <l_logical *>select_cpy.data
else:
MM = N
select_ptr = NULL
- if ((left and right and Q is not None and Z is not None) or
- (left and not right and Q is not None) or
- (right and not left and Z is not None)):
+ if Q is not None:
howmny = "B"
- backtr = True
else:
howmny = "A"
+ cdef np.ndarray[scalar, ndim=2] vl_r = None
+ cdef scalar *vl_r_ptr
if left:
- if backtr:
- vl_r = Q
+ if Q is not None and select is None:
+ vl_r = np.asfortranarray(Q.copy())
else:
- vl_r = np.empty((N, MM), dtype = np.float32, order='F')
- vl_r_ptr = <float *>vl_r.data
+ vl_r = np.empty((N, MM), dtype=T.dtype, order='F')
+ vl_r_ptr = <scalar *>vl_r.data
else:
vl_r_ptr = NULL
+ cdef np.ndarray[scalar, ndim=2] vr_r = None
+ cdef scalar *vr_r_ptr
if right:
- if backtr:
- vr_r = Z
+ if Q is not None and select is None:
+ vr_r = np.asfortranarray(Q.copy())
else:
- vr_r = np.empty((N, MM), dtype = np.float32, order='F')
- vr_r_ptr = <float *>vr_r.data
+ vr_r = np.empty((N, MM), dtype=T.dtype, order='F')
+ vr_r_ptr = <scalar *>vr_r.data
else:
vr_r_ptr = NULL
- f_lapack.stgevc_(side, howmny, select_ptr,
- &N, <float *>S.data, &N,
- <float *>T.data, &N,
- vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
- <float *>work.data, &info)
-
- assert info == 0, "Argument error in stgevc"
- assert MM == M, "Unexpected number of eigenvectors returned in stgevc"
+ # The actual calculation
+
+ if scalar is float:
+ lapack.strevc(side, howmny, select_ptr,
+ &N, <float *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <float *>work.data, &info)
+ elif scalar is double:
+ lapack.dtrevc(side, howmny, select_ptr,
+ &N, <double *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <double *>work.data, &info)
+ elif scalar is float_complex:
+ lapack.ctrevc(side, howmny, select_ptr,
+ &N, <float complex *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <float complex *>work.data, <float *>rwork.data, &info)
+ elif scalar is double_complex:
+ lapack.ztrevc(side, howmny, select_ptr,
+ &N, <double complex *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <double complex *>work.data, <double *>rwork.data, &info)
+
+ assert info == 0, "Argument error in trevc"
+ assert MM == M, "Unexpected number of eigenvectors returned in strevc"
- if not backtr:
+ if select is not None and Q is not None:
if left:
vl_r = np.asfortranarray(np.dot(Q, vl_r))
if right:
- vr_r = np.asfortranarray(np.dot(Z, vr_r))
+ vr_r = np.asfortranarray(np.dot(Q, vr_r))
- # If there are complex eigenvalues, we need to postprocess the eigenvectors
- if np.diagonal(S, -1).nonzero()[0].size:
- if left:
- vl = txevc_postprocess(np.complex64, S, vl_r, select)
- if right:
- vr = txevc_postprocess(np.complex64, S, vr_r, select)
- else:
- if left:
- vl = vl_r
- if right:
- vr = vr_r
+ cdef np.ndarray vl, vr
+ if left:
+ vl = vl_r
+ if right:
+ vr = vr_r
+ if scalar in floating:
+ # If there are complex eigenvalues, we need to postprocess the
+ # eigenvectors.
+ if scalar is float:
+ dtype = np.complex64
+ else:
+ dtype = np.complex128
+ if np.diagonal(T, -1).nonzero()[0].size:
+ if left:
+ vl = txevc_postprocess(dtype, T, vl_r, select)
+ if right:
+ vr = txevc_postprocess(dtype, T, vr_r, select)
if left and right:
return (vl, vr)
@@ -2170,229 +799,340 @@ def stgevc(np.ndarray[np.float32_t, ndim=2] S,
return vr
-def dtgevc(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,
- 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
+def gges(np.ndarray[scalar, ndim=2] A,
+ np.ndarray[scalar, ndim=2] B,
+ calc_q=True, calc_z=True, calc_ev=True):
+ cdef l_int N, sdim, info
- assert_fortran_mat(S, T, Q, Z)
+ # Check parameters
- N = S.shape[0]
- work = np.empty(6*N, dtype = np.float64)
+ assert_fortran_mat(A, B)
- if left and right:
- side = "B"
- elif left:
- side = "L"
- elif right:
- side = "R"
- else:
- return
+ if A.shape[0] != B.shape[1]:
+ raise ValueError("Expect square matrix A")
- backtr = False
+ if A.shape[0] != B.shape[0] or A.shape[0] != B.shape[1]:
+ raise ValueError("Shape of B is incompatible with matrix A")
- if select is not None:
- howmny = "S"
- MM = select.nonzero()[0].size
- # Correct for possible additional storage if a single complex
- # eigenvalue is selected.
- # For that: Figure out the positions of the 2x2 blocks.
- cmplxindx = np.diagonal(S, -1).nonzero()[0]
- for i in cmplxindx:
- if bool(select[i]) != bool(select[i+1]):
- MM += 1
+ # Allocate workspaces
- # select is overwritten in dtgevc
- select_cpy = np.array(select, dtype = f_lapack.l_logical_dtype,
- order = 'F')
- select_ptr = <l_logical *>select_cpy.data
+ N = A.shape[0]
+
+ cdef np.ndarray[scalar] alphar, alphai
+ if scalar in cmplx:
+ alphar = np.empty(N, dtype=A.dtype)
+ alphai = None
else:
- MM = N
- select_ptr = NULL
- if ((left and right and Q is not None and Z is not None) or
- (left and not right and Q is not None) or
- (right and not left and Z is not None)):
- howmny = "B"
- backtr = True
- else:
- howmny = "A"
+ alphar = np.empty(N, dtype=A.dtype)
+ alphai = np.empty(N, dtype=A.dtype)
- if left:
- if backtr:
- vl_r = Q
- else:
- vl_r = np.empty((N, MM), dtype = np.float64, order='F')
- vl_r_ptr = <double *>vl_r.data
+ 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 char *jobvsl
+ cdef scalar *vsl_ptr
+ cdef np.ndarray[scalar, ndim=2] vsl
+ if calc_q:
+ vsl = np.empty((N,N), dtype=A.dtype, order='F')
+ vsl_ptr = <scalar *>vsl.data
+ jobvsl = "V"
else:
- vl_r_ptr = NULL
+ vsl = None
+ vsl_ptr = NULL
+ jobvsl = "N"
- if right:
- if backtr:
- vr_r = Z
- else:
- vr_r = np.empty((N, MM), dtype = np.float64, order='F')
- vr_r_ptr = <double *>vr_r.data
+ cdef char *jobvsr
+ cdef scalar *vsr_ptr
+ cdef np.ndarray[scalar, ndim=2] vsr
+ if calc_z:
+ vsr = np.empty((N,N), dtype=A.dtype, order='F')
+ vsr_ptr = <scalar *>vsr.data
+ jobvsr = "V"
else:
- vr_r_ptr = NULL
+ vsr = None
+ vsr_ptr = NULL
+ jobvsr = "N"
+
+ # Workspace query
+ # Xgges expects &qwork as a <scalar *> (even though it's an integer)
+ cdef l_int lwork = -1
+ cdef scalar qwork
+
+ if scalar is float:
+ lapack.sgges(jobvsl, jobvsr, "N", NULL,
+ &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)
+ elif scalar is double:
+ lapack.dgges(jobvsl, jobvsr, "N", NULL,
+ &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)
+ elif scalar is float_complex:
+ lapack.cgges(jobvsl, jobvsr, "N", NULL,
+ &N, <float complex *>A.data, &N,
+ <float complex *>B.data, &N, &sdim,
+ <float complex *>alphar.data, <float complex *>beta.data,
+ vsl_ptr, &N, vsr_ptr, &N,
+ &qwork, &lwork, <float *>rwork.data, NULL, &info)
+ elif scalar is double_complex:
+ lapack.zgges(jobvsl, jobvsr, "N", NULL,
+ &N, <double complex *>A.data, &N,
+ <double complex *>B.data, &N, &sdim,
+ <double complex *>alphar.data, <double complex *>beta.data,
+ vsl_ptr, &N, vsr_ptr, &N,
+ &qwork, &lwork, <double *>rwork.data, NULL, &info)
+
+ assert info == 0, "Argument error in gges"
+
+ lwork = lwork_from_qwork(qwork)
+ cdef np.ndarray[scalar] work = np.empty(lwork, dtype=A.dtype)
+
+ # The actual calculation
+
+ if scalar is float:
+ lapack.sgges(jobvsl, jobvsr, "N", NULL,
+ &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)
+ elif scalar is double:
+ lapack.dgges(jobvsl, jobvsr, "N", NULL,
+ &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)
+ elif scalar is float_complex:
+ lapack.cgges(jobvsl, jobvsr, "N", NULL,
+ &N, <float complex *>A.data, &N,
+ <float complex *>B.data, &N, &sdim,
+ <float complex *>alphar.data, <float complex *>beta.data,
+ vsl_ptr, &N, vsr_ptr, &N,
+ <float complex *>work.data, &lwork,
+ <float *>rwork.data, NULL, &info)
+ elif scalar is double_complex:
+ lapack.zgges(jobvsl, jobvsr, "N", NULL,
+ &N, <double complex *>A.data, &N,
+ <double complex *>B.data, &N, &sdim,
+ <double complex *>alphar.data, <double complex *>beta.data,
+ vsl_ptr, &N, vsr_ptr, &N,
+ <double complex *>work.data, &lwork,
+ <double *>rwork.data, NULL, &info)
- f_lapack.dtgevc_(side, howmny, select_ptr,
- &N, <double *>S.data, &N,
- <double *>T.data, &N,
- vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
- <double *>work.data, &info)
+ if info > 0:
+ raise LinAlgError("QZ iteration failed to converge in gges")
- assert info == 0, "Argument error in dtgevc"
- assert MM == M, "Unexpected number of eigenvectors returned in dtgevc"
+ assert info == 0, "Argument error in gges"
- if not backtr:
- if left:
- vl_r = np.asfortranarray(np.dot(Q, vl_r))
- if right:
- vr_r = np.asfortranarray(np.dot(Z, vr_r))
+ # Real inputs possibly produce complex output
+ cdef np.ndarray alpha = maybe_complex[scalar](0, alphar, alphai)
- # If there are complex eigenvalues, we need to postprocess the
- # eigenvectors.
- if np.diagonal(S, -1).nonzero()[0].size:
- if left:
- vl = txevc_postprocess(np.complex128, S, vl_r, select)
- if right:
- vr = txevc_postprocess(np.complex128, S, vr_r, select)
- else:
- if left:
- vl = vl_r
- if right:
- vr = vr_r
+ return filter_args((True, True, calc_q, calc_z, calc_ev, calc_ev),
+ (A, B, vsl, vsr, alpha, beta))
- if left and right:
- return (vl, vr)
- elif left:
- return vl
- else:
- return vr
+def tgsen(np.ndarray[l_logical] select,
+ np.ndarray[scalar, ndim=2] S,
+ np.ndarray[scalar, ndim=2] T,
+ np.ndarray[scalar, ndim=2] Q,
+ np.ndarray[scalar, ndim=2] Z,
+ calc_ev=True):
+ cdef l_int ijob = 0
+ cdef l_int N, M, lwork, liwork, info
-def ctgevc(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,
- 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
+ # Check parameters
+
+ if ((S.shape[0] != S.shape[1] or T.shape[0] != T.shape[1] or
+ S.shape[0] != T.shape[0]) or
+ (Q is not None and (Q.shape[0] != Q.shape[1] or
+ S.shape[0] != Q.shape[0])) or
+ (Z is not None and (Z.shape[0] != Z.shape[1] or
+ S.shape[0] != Z.shape[0]))):
+ raise ValueError("Invalid Schur decomposition as input")
assert_fortran_mat(S, T, Q, Z)
+ # Allocate workspaces
+
N = S.shape[0]
- work = np.empty(2*N, dtype = np.complex64)
- rwork = np.empty(2*N, dtype = np.float32)
- if left and right:
- side = "B"
- elif left:
- side = "L"
- elif right:
- side = "R"
+ cdef np.ndarray[scalar] alphar, alphai
+ if scalar in cmplx:
+ alphar = np.empty(N, dtype=S.dtype)
+ alphai = None
else:
- return
+ alphar = np.empty(N, dtype=S.dtype)
+ alphai = np.empty(N, dtype=S.dtype)
- backtr = False
+ cdef np.ndarray[scalar] beta
+ beta = np.empty(N, dtype=S.dtype)
- if select is not None:
- howmny = "S"
- MM = select.nonzero()[0].size
- select_ptr = <l_logical *>select.data
+ cdef l_logical wantq
+ cdef scalar *q_ptr
+ if Q is not None:
+ wantq = 1
+ q_ptr = <scalar *>Q.data
else:
- MM = N
- select_ptr = NULL
- if ((left and right and Q is not None and Z is not None) or
- (left and not right and Q is not None) or
- (right and not left and Z is not None)):
- howmny = "B"
- backtr = True
- else:
- howmny = "A"
+ wantq = 0
+ q_ptr = NULL
- if left:
- if backtr:
- vl = Q
- else:
- vl = np.empty((N, MM), dtype = np.complex64, order='F')
- vl_ptr = <float complex *>vl.data
+ cdef l_logical wantz
+ cdef scalar *z_ptr
+ if Z is not None:
+ wantz = 1
+ z_ptr = <scalar *>Z.data
else:
- vl_ptr = NULL
+ wantz = 0
+ z_ptr = NULL
- if right:
- if backtr:
- vr = Z
- else:
- vr = np.empty((N, MM), dtype = np.complex64, order='F')
- vr_ptr = <float complex *>vr.data
- else:
- vr_ptr = NULL
+ # Workspace query
+ # Xtgsen expects &qwork as a <scalar *> (even though it's an integer)
+ lwork = -1
+ liwork = -1
+ cdef scalar qwork
+ cdef l_int qiwork
+
+ if scalar is float:
+ 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)
+ elif scalar is double:
+ lapack.dtgsen(&ijob, &wantq, &wantz, <l_logical *>select.data,
+ &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)
+ elif scalar is float_complex:
+ lapack.ctgsen(&ijob, &wantq, &wantz, <l_logical *>select.data,
+ &N, <float complex *>S.data, &N,
+ <float complex *>T.data, &N,
+ <float complex *>alphar.data, <float complex *>beta.data,
+ q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL,
+ &qwork, &lwork, &qiwork, &liwork, &info)
+ elif scalar is double_complex:
+ lapack.ztgsen(&ijob, &wantq, &wantz, <l_logical *>select.data,
+ &N, <double complex *>S.data, &N,
+ <double complex *>T.data, &N,
+ <double complex *>alphar.data, <double complex *>beta.data,
+ q_ptr, &N, z_ptr, &N, &M, NULL, NULL, NULL,
+ &qwork, &lwork, &qiwork, &liwork, &info)
+
+ assert info == 0, "Argument error in tgsen"
+
+ lwork = lwork_from_qwork(qwork)
+ cdef np.ndarray[scalar] work = np.empty(lwork, dtype=S.dtype)
- f_lapack.ctgevc_(side, howmny, select_ptr,
- &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)
+ liwork = qiwork
+ cdef np.ndarray[l_int] iwork = np.empty(liwork, dtype=int_dtype)
+
+ # The actual calculation
+
+ if scalar is float:
+ 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)
+ elif scalar is double:
+ lapack.dtgsen(&ijob, &wantq, &wantz, <l_logical *>select.data,
+ &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)
+ elif scalar is float_complex:
+ lapack.ctgsen(&ijob, &wantq, &wantz, <l_logical *>select.data,
+ &N, <float complex *>S.data, &N,
+ <float complex *>T.data, &N,
+ <float complex *>alphar.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)
+ elif scalar is double_complex:
+ lapack.ztgsen(&ijob, &wantq, &wantz, <l_logical *>select.data,
+ &N, <double complex *>S.data, &N,
+ <double complex *>T.data, &N,
+ <double complex *>alphar.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)
- assert info == 0, "Argument error in ctgevc"
- assert MM == M, "Unexpected number of eigenvectors returned in ctgevc"
+ if info > 0:
+ raise LinAlgError("Reordering failed; problem is very ill-conditioned")
- if not backtr:
- if left:
- vl = np.asfortranarray(np.dot(Q, vl))
- if right:
- vr = np.asfortranarray(np.dot(Z, vr))
+ assert info == 0, "Argument error in tgsen"
- if left and right:
- return (vl, vr)
- elif left:
- return vl
- else:
- return vr
+ # Real inputs possibly produce complex output
+ cdef np.ndarray alpha = maybe_complex[scalar](0, alphar, alphai)
+
+ return filter_args((True, True, Q is not None, Z is not None,
+ calc_ev, calc_ev),
+ (S, T, Q, Z, alpha, beta))
-def ztgevc(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,
- np.ndarray[l_logical] select=None,
- left=False, right=True):
+def tgevc(np.ndarray[scalar, ndim=2] S,
+ np.ndarray[scalar, ndim=2] T,
+ np.ndarray[scalar, ndim=2] Q,
+ np.ndarray[scalar, ndim=2] Z,
+ np.ndarray[l_logical] select,
+ 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
+
+ # Check parameters
+
+ if ((S.shape[0] != S.shape[1] or T.shape[0] != T.shape[1] or
+ S.shape[0] != T.shape[0]) or
+ (Q is not None and (Q.shape[0] != Q.shape[1] or
+ S.shape[0] != Q.shape[0])) or
+ (Z is not None and (Z.shape[0] != Z.shape[1] or
+ S.shape[0] != Z.shape[0]))):
+ raise ValueError("Invalid Schur decomposition as input")
assert_fortran_mat(S, T, Q, Z)
+ # Allocate workspaces
+
N = S.shape[0]
- work = np.empty(2*N, dtype = np.complex128)
- rwork = np.empty(2*N, dtype = np.float64)
+ cdef np.ndarray[scalar] work
+ if scalar in floating:
+ work = np.empty(6 * N, dtype=S.dtype)
+ else:
+ work = np.empty(2 * N, dtype=S.dtype)
+
+ cdef np.ndarray rwork = None
+ if scalar is float_complex:
+ rwork = np.empty(2 * N, dtype=np.float32)
+ elif scalar is double_complex:
+ rwork = np.empty(2 * N, dtype=np.float64)
+
+ cdef char *side
if left and right:
side = "B"
elif left:
@@ -2402,12 +1142,26 @@ def ztgevc(np.ndarray[np.complex128_t, ndim=2] S,
else:
return
- backtr = False
+ cdef l_logical backtr = False
+ cdef char *howmny
+ cdef np.ndarray[l_logical] select_cpy = None
+ cdef l_logical *select_ptr
if select is not None:
howmny = "S"
MM = select.nonzero()[0].size
- select_ptr = <l_logical *>select.data
+ # Correct for possible additional storage if a single complex
+ # eigenvalue is selected.
+ # For that: Figure out the positions of the 2x2 blocks.
+ cmplxindx = np.diagonal(S, -1).nonzero()[0]
+ for i in cmplxindx:
+ if bool(select[i]) != bool(select[i+1]):
+ MM += 1
+
+ # select is overwritten in tgevc
+ select_cpy = np.array(select, dtype=logical_dtype,
+ order = 'F')
+ select_ptr = <l_logical *>select_cpy.data
else:
MM = N
select_ptr = NULL
@@ -2419,38 +1173,78 @@ def ztgevc(np.ndarray[np.complex128_t, ndim=2] S,
else:
howmny = "A"
+ cdef np.ndarray[scalar, ndim=2] vl_r
+ cdef scalar *vl_r_ptr
if left:
if backtr:
- vl = Q
+ vl_r = Q
else:
- vl = np.empty((N, MM), dtype = np.complex128, order='F')
- vl_ptr = <double complex *>vl.data
+ vl_r = np.empty((N, MM), dtype=S.dtype, order='F')
+ vl_r_ptr = <scalar *>vl_r.data
else:
- vl_ptr = NULL
+ vl_r_ptr = NULL
+ cdef np.ndarray[scalar, ndim=2] vr_r
+ cdef scalar *vr_r_ptr
if right:
if backtr:
- vr = Z
+ vr_r = Z
else:
- vr = np.empty((N, MM), dtype = np.complex128, order='F')
- vr_ptr = <double complex *>vr.data
+ vr_r = np.empty((N, MM), dtype=S.dtype, order='F')
+ vr_r_ptr = <scalar *>vr_r.data
else:
- vr_ptr = NULL
-
- f_lapack.ztgevc_(side, howmny, select_ptr,
- &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)
+ vr_r_ptr = NULL
- assert info == 0, "Argument error in ztgevc"
- assert MM == M, "Unexpected number of eigenvectors returned in ztgevc"
+ if scalar is float:
+ lapack.stgevc(side, howmny, select_ptr,
+ &N, <float *>S.data, &N,
+ <float *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <float *>work.data, &info)
+ elif scalar is double:
+ lapack.dtgevc(side, howmny, select_ptr,
+ &N, <double *>S.data, &N,
+ <double *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <double *>work.data, &info)
+ elif scalar is float_complex:
+ lapack.ctgevc(side, howmny, select_ptr,
+ &N, <float complex *>S.data, &N,
+ <float complex *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <float complex *>work.data, <float *>rwork.data, &info)
+ elif scalar is double_complex:
+ lapack.ztgevc(side, howmny, select_ptr,
+ &N, <double complex *>S.data, &N,
+ <double complex *>T.data, &N,
+ vl_r_ptr, &N, vr_r_ptr, &N, &MM, &M,
+ <double complex *>work.data, <double *>rwork.data, &info)
+
+ assert info == 0, "Argument error in tgevc"
+ assert MM == M, "Unexpected number of eigenvectors returned in tgevc"
if not backtr:
if left:
- vl = np.asfortranarray(np.dot(Q, vl))
+ vl_r = np.asfortranarray(np.dot(Q, vl_r))
if right:
- vr = np.asfortranarray(np.dot(Z, vr))
+ vr_r = np.asfortranarray(np.dot(Z, vr_r))
+
+ # If there are complex eigenvalues, we need to postprocess the eigenvectors
+ cdef np.ndarray vl, vr
+ if left:
+ vl = vl_r
+ if right:
+ vr = vr_r
+ if scalar in floating:
+ if scalar is float:
+ dtype = np.complex64
+ else:
+ dtype = np.complex128
+ if np.diagonal(S, -1).nonzero()[0].size:
+ if left:
+ vl = txevc_postprocess(dtype, S, vl_r, select)
+ if right:
+ vr = txevc_postprocess(dtype, S, vr_r, select)
if left and right:
return (vl, vr)
@@ -2480,8 +1274,7 @@ def prepare_for_lapack(overwrite, *args):
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.
+ Returns a list of properly converted arrays.
"""
# Make sure we have NumPy arrays
@@ -2502,18 +1295,10 @@ def prepare_for_lapack(overwrite, *args):
# kind.
dtype = np.common_type(*[arr for arr, ovwrt in mats if arr is not None])
- if dtype == np.float32:
- lapacktype = 's'
- elif dtype == np.float64:
- lapacktype = 'd'
- elif dtype == np.complex64:
- lapacktype = 'c'
- elif dtype == np.complex128:
- lapacktype = 'z'
- else:
+ if dtype not in (np.float32, np.float64, np.complex64, np.complex128):
raise AssertionError("Unexpected data type from common_type")
- ret = [ lapacktype ]
+ ret = []
for npmat, ovwrt in mats:
# Now make sure that the array is contiguous, and copy if necessary.
if npmat is not None:
@@ -2538,4 +1323,7 @@ def prepare_for_lapack(overwrite, *args):
ret.append(npmat)
- return tuple(ret)
+ if len(ret) == 1:
+ return ret[0]
+ else:
+ return tuple(ret)
diff --git a/setup.py b/setup.py
index 197505db81f8a782506519d5a4f70687f713aa84..04d92ea7a23deb14abd64ef07aa4cc7cdb7ddfd4 100755
--- a/setup.py
+++ b/setup.py
@@ -372,8 +372,9 @@ def search_mumps():
# Conda (via conda-forge).
# TODO: remove dependency libs (scotch, metis...) when conda-forge
# packaged mumps/scotch are built as properly linked shared libs
+ # 'openblas' provides Lapack and BLAS symbols
['zmumps', 'mumps_common', 'metis', 'esmumps', 'scotch',
- 'scotcherr', 'mpiseq'],
+ 'scotcherr', 'mpiseq', 'openblas'],
]
common_libs = ['pord', 'gfortran']
@@ -384,34 +385,7 @@ def search_mumps():
return []
-def search_lapack():
- """Return the BLAS variant that is installed."""
- lib_sets = [
- # Debian
- ['blas', 'lapack'],
- # Conda (via conda-forge). Openblas contains lapack symbols
- ['openblas', 'gfortran'],
- ]
-
- for libs in lib_sets:
- found_libs = search_libs(libs)
- if found_libs:
- return found_libs
-
- print('Error: BLAS/LAPACK are required but were not found.',
- file=sys.stderr)
- sys.exit(1)
-
-
def configure_special_extensions(exts, build_summary):
- #### Special config for LAPACK.
- lapack = exts['kwant.linalg.lapack']
- if 'libraries' in lapack:
- build_summary.append('User-configured LAPACK and BLAS')
- else:
- lapack['libraries'] = search_lapack()
- build_summary.append('Default LAPACK and BLAS')
-
#### Special config for MUMPS.
mumps = exts['kwant.linalg._mumps']
if 'libraries' in mumps:
@@ -426,12 +400,6 @@ def configure_special_extensions(exts, build_summary):
del exts['kwant.linalg._mumps']
build_summary.append('No MUMPS support')
- if mumps:
- # Copy config from LAPACK.
- for key, value in lapack.items():
- if key not in ['sources', 'depends']:
- mumps.setdefault(key, []).extend(value)
-
return exts
@@ -550,8 +518,7 @@ def main():
'kwant/graph/c_slicer/partitioner.h',
'kwant/graph/c_slicer/slicer.h'])),
('kwant.linalg.lapack',
- dict(sources=['kwant/linalg/lapack.pyx'],
- depends=['kwant/linalg/f_lapack.pxd'])),
+ dict(sources=['kwant/linalg/lapack.pyx'])),
('kwant.linalg._mumps',
dict(sources=['kwant/linalg/_mumps.pyx'],
depends=['kwant/linalg/cmumps.pxd']))])
@@ -567,8 +534,7 @@ def main():
for ext in exts.values():
ext.setdefault('include_dirs', []).append(numpy_include)
- aliases = [('lapack', 'kwant.linalg.lapack'),
- ('mumps', 'kwant.linalg._mumps')]
+ aliases = [('mumps', 'kwant.linalg._mumps')]
init_cython()