In this chapter we briefly describe the Python calling sequences. For
further details on the underlying LAPACK functions we refer to the LAPACK
Users’ Guide and manual pages.

The BLAS conventional storage scheme of the section Matrix Classes
is used. As in the previous chapter, we omit from the function definitions
less important arguments that are useful for selecting submatrices. The
complete definitions are documented in the docstrings in the source code.

The arguments A and B must have the same type ('d'
or 'z'). On entry, B contains the right-hand side
; on exit it contains the solution . The optional
argument ipiv is an integer matrix of length at least .
If ipiv is provided, then gesv solves the system, replaces
A with the triangular factors in an LU factorization, and returns
the permutation matrix in ipiv. If ipiv is not specified,
then gesv solves the system but does not return the LU
factorization and does not modify A.

The argument ipiv is an integer matrix of length at least
min{, }. On exit, the lower triangular part of
A is replaced by , the upper triangular part by ,
and the permutation matrix is returned in ipiv.

where and are real or complex matrices, with
by and banded with
subdiagonals.

The arguments A and B must have the same type ('d'
or 'z'). On entry, B contains the right-hand side
; on exit it contains the solution . The optional
argument ipiv is an integer matrix of length at least .
If ipiv is provided, then A must have
rows. On entry the diagonals of are stored in rows
to of A, using the BLAS
format for general band matrices (see the section
Matrix Classes). On exit, the factorization is returned in
A and ipiv. If ipiv is not provided, then A must have
rows. On entry the diagonals of are
stored in the rows of A, following the standard BLAS format for
general band matrices. In this case, gbsv does not modify
A and does not return the factorization.

LU factorization of a general by real or complex
band matrix with subdiagonals.

The matrix is stored using the BLAS format for general band matrices
(see the section Matrix Classes), by providing the diagonals
(stored as rows of a by matrix A),
the number of rows , and the number of subdiagonals
. The argument ipiv is an integer matrix of length at
least min{, }. On exit, A and ipiv contain
the details of the factorization.

The subdiagonal of is stored as a matrix dl of length
, the diagonal is stored as a matrix d of length
, and the superdiagonal is stored as a matrix du of
length . The four arguments must have the same type
('d' or 'z'). On exit dl, d, du are
overwritten with the details of the LU factorization of .
On entry, B contains the right-hand side ; on exit it
contains the solution .

The subdiagonal of is stored as a matrix dl of length
, the diagonal is stored as a matrix d of length
, and the superdiagonal is stored as a matrix du of length
. dl, d and du must have the same type.
du2 is a matrix of length , and of the same type as
dl. ipiv is an 'i' matrix of length .
On exit, the five arguments contain the details of the factorization.

The arguments dl, d, du, du2, and ipiv contain
the details of the LU factorization as returned by
gttrf.
On entry, B contains the right-hand side ; on exit it
contains the solution . B must have the same type as
the other arguments.

where is a real symmetric or complex Hermitian positive
definite band matrix.

On entry, the diagonals of are stored in A, using the
BLAS format for symmetric or Hermitian band matrices (see
section Matrix Classes). On exit, B is replaced by the
solution, and A is overwritten with the Cholesky factor (in the
BLAS format for triangular band matrices). The matrices A and
B must have the same type ('d' or 'z').

where is an by positive definite real
symmetric or complex Hermitian tridiagonal matrix.

The diagonal of is stored as a 'd' matrix d of
length and its subdiagonal as a 'd' or 'z'
matrix e of length . The arguments e and B
must have the same type. On exit d contains the diagonal elements
of in the
LDLT
or
LDLH
factorization of , and
e contains the subdiagonal elements of the unit lower bidiagonal
matrix . B is overwritten with the solution .
Raises an ArithmeticError if the matrix is singular.

On entry, the argument d is a 'd' matrix with the diagonal
elements of . The argument e is 'd' or
'z' matrix containing the subdiagonal of . On exit
d contains the diagonal elements of , and e contains
the subdiagonal elements of the unit lower bidiagonal matrix .

where is an by positive definite real
symmetric or complex Hermitian tridiagonal matrix, given its
LDLT
or
LDLH
factorization.

The argument d is the diagonal of the diagonal matrix .
The argument uplo only matters for complex matrices. If uplo
is 'L', then on exit e contains the subdiagonal elements
of the unit bidiagonal matrix . If uplo is 'U',
then e contains the complex conjugates of the elements of the unit
bidiagonal matrix . On exit, B is overwritten with the
solution . B must have the same type as e.

On exit, B is replaced by the solution. The matrices A and
B must have the same type ('d' or 'z'). The
optional argument ipiv is an integer matrix of length at least
equal to . If ipiv is provided, sysv solves the
system and returns the factorization in A and ipiv. If
ipiv is not specified, sysv solves the system but does not
return the factorization and does not modify A.

On exit, B is replaced by the solution. The matrices A and
B must have the same type ('d' or 'z'). The
optional argument ipiv is an integer matrix of length at least
. If ipiv is provided, then hesv solves the
system and returns the factorization in A and ipiv. If
ipiv is not specified, then hesv solves the system but does
not return the factorization and does not modify A.

A and B must have the same typecode ('d' or
'z'). trans = 'T' is not allowed if A is
complex. On exit, the solution is stored as the leading
submatrix of B. The matrix A is overwritten with details of
the QR or the LQ factorization of .

If is by , then is by
and orthogonal/unitary, and is by
and upper triangular (if is greater than or equal
to ), or upper trapezoidal (if is less than or
equal to ).

tau is a matrix of the same type as A and of length
min{, }. On exit, is stored in the upper
triangular/trapezoidal part of A. The matrix is stored
as a product of min{, } elementary reflectors in
the first min{, } columns of A and in tau.

If is by , then is by
and orthogonal/unitary, and is by
and lower triangular (if is less than or equal to
), or lower trapezoidal (if is greater than or equal
to ).

tau is a matrix of the same type as A and of length
min{, }. On exit, is stored in the lower
triangular/trapezoidal part of A. The matrix is stored
as a product of min{, } elementary reflectors in the
first min{, } rows of A and in tau.

If is by , then is
by and orthogonal/unitary, and is by
and upper triangular (if is greater than or equal
to ), or upper trapezoidal (if is less than or equal
to ).

tau is a matrix of the same type as A and of length
min{, }. jpvt is an integer matrix of
length . On entry, if jpvt[k] is nonzero, then
column of is permuted to the front of .
Otherwise, column is a free column.

On exit, jpvt contains the permutation : the operation
is equivalent to A[:,jpvt-1]. is stored
in the upper triangular/trapezoidal part of A. The matrix
is stored as a product of min{, }
elementary reflectors in the first min{,:math:n} columns
of A and in tau.

In most applications, the matrix is not needed explicitly, and
it is sufficient to be able to make products with or its
transpose. The functions unmqr and
ormqr multiply a matrix
with the orthogonal matrix computed by
geqrf.

If A is by , then is square of order
and orthogonal or unitary. is stored in the first
min{, } columns of A and in tau as a
product of min{, } elementary reflectors, as
computed by geqrf.
The matrices A, tau, and C
must have the same type. trans = 'T' is only allowed if
the typecode is 'd'.

If A is by , then is square of order
and orthogonal or unitary. is stored in the first
min{, } rows of A and in tau as a product of
min{, } elementary reflectors, as computed by
gelqf.
The matrices A, tau, and C must have the
same type. trans = 'T' is only allowed if the typecode
is 'd'.

If A has size by , and tau has length
, then, on entry, the first k columns of the matrix A
and the entries of tau contai an unitary or orthogonal matrix
of order , as computed by
geqrf. On exit,
the first min{, } columns of are contained
in the leading columns of A.

If A has size by , and tau has length
, then, on entry, the first k rows of the matrix A
and the entries of tau contain a unitary or orthogonal matrix
of order , as computed by
gelqf.
On exit, the first min{, } rows of are
contained in the leading rows of A.

W is a real matrix of length at least . On exit, W
contains the eigenvalues in ascending order. If jobz is
'V', the eigenvectors are also computed and returned in A.
If jobz is 'N', the eigenvectors are not returned and the
contents of A are destroyed.

Computes selected eigenvalues and eigenvectors of a real symmetric
matrix of order .

W is a real matrix of length at least . On exit, W
contains the eigenvalues in ascending order. If range is
'A', all the eigenvalues are computed. If range is
'I', eigenvalues through are
computed, where . If range is
'V', the eigenvalues in the interval are
computed.

If jobz is 'V', the (normalized) eigenvectors are
computed, and returned in Z. If jobz is 'N', the
eigenvectors are not computed. In both cases, the contents of A
are destroyed on exit.

Z is optional (and not referenced) if jobz is 'N'.
It is required if jobz is 'V' and must have at least
columns if range is 'A' or 'V' and at
least columns if range is 'I'.

Solves the generalized eigenproblem (2) for real symmetric
matrices of order , stored in real matrices A and B.
itype is an integer with possible values 1, 2, 3, and specifies
the type of eigenproblem. W is a real matrix of length at least
. On exit, it contains the eigenvalues in ascending order.
On exit, B contains the Cholesky factor of . If jobz
is 'V', the eigenvectors are computed and returned in A.
If jobz is 'N', the eigenvectors are not returned and the
contents of A are destroyed.

S is a real matrix of length at least min{, }.
On exit, its first min{, } elements are the
singular values in descending order.

The argument jobu controls how many left singular vectors are
computed. The possible values are 'N', 'A',
'S' and 'O'. If jobu is 'N', no left
singular vectors are computed. If jobu is 'A', all left
singular vectors are computed and returned as columns of U.
If jobu is 'S', the first min{, } left
singular vectors are computed and returned as columns of U.
If jobu is 'O', the first min{, } left
singular vectors are computed and returned as columns of A.
The argument U is None(if jobu is 'N'
or 'A') or a matrix of the same type as A.

The argument jobvt controls how many right singular vectors are
computed. The possible values are 'N', 'A',
'S' and 'O'. If jobvt is 'N', no right
singular vectors are computed. If jobvt is 'A', all right
singular vectors are computed and returned as rows of Vt.
If jobvt is 'S', the first min{, }
right singular vectors are computed and their (conjugate) transposes
are returned as rows of Vt. If jobvt is 'O', the
first min{, } right singular vectors are computed
and their (conjugate) transposes are returned as rows of A.
Note that the (conjugate) transposes of the right singular vectors
(i.e., the matrix ) are returned in Vt or A.
The argument Vt can be None (if jobvt is 'N'
or 'A') or a matrix of the same type as A.

Singular value decomposition of a real or complex by
matrix.. This function is based on a divide-and-conquer
algorithm and is faster than gesvd.

S is a real matrix of length at least min{, }.
On exit, its first min{, } elements are the
singular values in descending order.

The argument jobz controls how many singular vectors are computed.
The possible values are 'N', 'A', 'S' and
'O'. If jobz is 'N', no singular vectors are
computed. If jobz is 'A', all left singular
vectors are computed and returned as columns of U and all
right singular vectors are computed and returned as rows of
Vt. If jobz is 'S', the first
min{, } left and right singular vectors are computed
and returned as columns of U and rows of Vt.
If jobz is 'O' and is greater than or equal
to , the first left singular vectors are returned as
columns of A and the right singular vectors are returned
as rows of Vt. If jobz is 'O' and is less
than , the left singular vectors are returned as
columns of U and the first right singular vectors are
returned as rows of A. Note that the (conjugate) transposes of
the right singular vectors are returned in Vt or A.

The argument U can be None (if jobz is 'N'
or 'A' of jobz is 'O' and is greater
than or equal to ) or a matrix of the same type as A.
The argument Vt can be None(if jobz is 'N'
or 'A' or jobz is 'O' and :math`m` is less than
) or a matrix of the same type as A.

If is real, the matrix of Schur vectors is
orthogonal, and is a real upper quasi-triangular matrix with
1 by 1 or 2 by 2 diagonal blocks. The 2 by 2 blocks correspond to
complex conjugate pairs of eigenvalues of .
If is complex, the matrix of Schur vectors is
unitary, and is a complex upper triangular matrix with the
eigenvalues of on the diagonal.

The optional argument w is a complex matrix of length at least
. If it is provided, the eigenvalues of A are returned
in w. The optional argument V is an by
matrix of the same type as A. If it is provided, then the Schur
vectors are returned in V.

The argument select is an optional ordering routine. It must be a
Python function that can be called as f(s) with a complex
argument s, and returns True or False. The
eigenvalues for which select returns True will be selected
to appear first along the diagonal. (In the real Schur factorization,
if either one of a complex conjugate pair of eigenvalues is selected,
then both are selected.)

On exit, A is replaced with the matrix . The function
gees returns an integer equal to the number of eigenvalues
that were selected by the ordering routine. If select is
None, then gees returns 0.

If and are real, then the matrices of left and
right Schur vectors and are orthogonal,
is a real upper quasi-triangular matrix with 1 by 1 or 2 by
2 diagonal blocks, and is a real triangular matrix with
nonnegative diagonal. The 2 by 2 blocks along the diagonal of
correspond to complex conjugate pairs of generalized
eigenvalues of , . If and are
complex, the matrices of left and right Schur vectors and
are unitary, is complex upper triangular, and
is complex upper triangular with nonnegative real diagonal.

The optional arguments a and b are 'z' and
'd' matrices of length at least . If these are
provided, the generalized eigenvalues of A, B are returned in
a and b. (The generalized eigenvalues are the ratios
a[k]/b[k].) The optional arguments Vl and Vr are
by matrices of the same type as A and B.
If they are provided, then the left Schur vectors are returned in
Vl and the right Schur vectors are returned in Vr.

The argument select is an optional ordering routine. It must be
a Python function that can be called as f(x,y) with a complex
argument x and a real argument y, and returns True or
False. The eigenvalues for which select returns
True will be selected to appear first on the diagonal.
(In the real Schur factorization, if either one of a complex conjugate
pair of eigenvalues is selected, then both are selected.)

On exit, A is replaced with the matrix and B is
replaced with the matrix . The function gges returns
an integer equal to the number of eigenvalues that were selected by
the ordering routine. If select is None, then
gges returns 0.

As an example, we compute the generalized complex Schur form of the
matrix of the previous example, and

In the code below we solve the problem using Newton’s method. At each
iteration the Newton direction is computed by solving a positive definite
set of linear equations

(where has rows ), and a suitable step size is
determined by a backtracking line search.

We use the level-3 BLAS function blas.syrk to
form the Hessian
matrix and the LAPACK function posv to
solve the Newton system.
The code can be further optimized by replacing the matrix-vector products
with the level-2 BLAS function blas.gemv.