SF_LINALG

Description

SciFortran module for linear algebra

Quick access

Used modules

sf_blacs: SciFortran BLACS module (used by SF_LINALG)

Subroutines and functions

interface sf_linalg/eig(a, eval, evec[, jobvl, jobvr])

This subroutine computes the eigenvalues and the right eigenvectors of a general (non-symmetric real or non-hermitian complex) square matrix A, \(A v_j = \lambda_j v_j\), using the LAPACK routines dgeev (real, deig) and zgeev (complex, zeig).

The input matrix is not modified. The eigenvalues Eval are complex and are returned in the order given by LAPACK. The eigenvectors are stored column-wise in Evec, Evec(:,j) being the eigenvector associated to Eval(j). The right eigenvectors are computed unless jobvr='N'. The left eigenvectors are not available: the program stops if jobvl='V', if A is not square, or if Evec does not have shape [n,n].

Parameters:

a (•, •) [real/complex, in] – real matrix [n,n], not modified
eval (•) [complex, out] – eigenvalues, size n
evec (•, •) [complex, out] – right eigenvectors, [n,n]; Evec(:,j) <-> Eval(j)

Options:

jobvl [character(len=1)] – optional, ‘N’ (default); ‘V’ is not supported
jobvr [character(len=1)] – optional, ‘V’ (default) computes the right eigenvectors, ‘N’ does not

interface sf_linalg/eigh(am, bm, lam, c, a, e, d, u[, method, jobz, uplo, vl, vu, il, iu, tol, ev, irange, vrange])

This subroutine solves the eigenvalue problem of real symmetric or complex hermitian matrices. The eigenvalues are real and are returned in ascending order. The specific procedures cover

the generalized problem \(A c = \lambda B c\), with Am symmetric/hermitian and Bm positive definite (dsygvd, zhegvd), returning the eigenvalues lam and the eigenvectors c, stored column-wise and normalized as \(c^H B c = 1\). Only the lower triangles of Am and Bm are referenced and the matrices are not modified
the simple problem \(A v = E v\), solved in place with the LAPACK driver chosen by method: dsyevd (default), dsyevr, dsyev or dsyevx, and the zhe... counterparts for hermitian matrices. If jobz='V' (default) the eigenvectors overwrite the columns of A, otherwise A is destroyed. The triangle uplo (default 'U') is used. A subset of the spectrum can be selected, with the r and x drivers only, by value (vl, vu) or by index (il, iu)
a real symmetric tridiagonal matrix, given by its main diagonal D and off-diagonal U and diagonalized with dstevr. D and U are overwritten, and D contains the eigenvalues on output. The eigenvectors are computed only if Ev is present, and a subset of the spectrum can be selected by index with Irange or by value with Vrange

The program stops if the shapes of the arguments are inconsistent, or if the LAPACK routine returns an error code.

Parameters:

am (•, •) [real/complex, in] – LHS matrix of Am c = lam Bm c, [n,n]
bm (•, •) [real/complex, in] – Bmt temporaries overwritten by dsygvd
lam (•) [real, out] – eigenvalues, ascending, size n
c (•, •) [real/complex, out] – eigenvectors, c(i,j) = i-th component of j-th vector, [n,n]
a (•, •) [real/complex, inout/in,required] – in: symmetric matrix; out: eigenvectors (columns) if jobz=V
e (size(a, 2)) [real, inout] – eigenvalues, ascending, size n (zero beyond the mE found)
d (•) [real, in,required] – in: main diagonal, size n; out: eigenvalues (first m valid)
u (max(1, -1 + size(d))) [real] – off-diagonal, size n-1; destroyed on exit

Options:

method [character(len=*)] – optional: ‘dsyevd’ (default), ‘dsyevr’, ‘dsyev’, ‘dsyevx’
jobz [character(len=1)] – optional: ‘V’ (default) eigenvectors too, ‘N’ eigenvalues only
uplo [character(len=1)] – optional: ‘U’ (default) or ‘L’, triangle of A to be used
vl [real] – optional: lower bound of the eigenvalue interval (vl,vu]
vu [real] – optional: upper bound of the eigenvalue interval (vl,vu]
il [integer] – optional: index of the first eigenvalue wanted
iu [integer] – optional: index of the last eigenvalue wanted
tol [real] – optional: absolute eigenvalue tolerance (dsyevr/dsyevx), default 0
ev (•, •) [real] – optional out: eigenvectors, [n,m]
irange (2) [integer] – optional: [il,iu] indices of the first and last eigenvalue wanted
vrange (2) [integer] – optional: [vl,vu] eigenvalue interval (vl,vu]

interface sf_linalg/p_eigh(a, w, nblock[, method, jobz, uplo, vl, vu, il, iu, tol])

This subroutine is the distributed-memory counterpart of eigh for the simple problem, for real symmetric (p_deigh_simple) and complex hermitian (p_zeigh_simple) matrices, using ScaLAPACK. It is available only when compiling with _SCALAPACK.

The matrix A, a full copy of which is expected on every process, is distributed on the BLACS process grid of the module in a block-cyclic layout with square blocks of size Nblock, and diagonalized with the driver chosen by method: PDSYEVR (default), PDSYEV, PDSYEVD or PDSYEVX. The eigenvalues are returned in ascending order in W, and if jobz='V' the eigenvectors are gathered back into A, column-wise.

Parameters:

a (•, •) [real/complex, inout/in,required] – in: symmetric matrix; out: eigenvectors (columns) if jobz=V
w (size(a, 2)) [real, inout] – eigenvalues, ascending, size n
nblock [integer] – block size of the block-cyclic distribution

Options:

method [character(len=*)] – optional: ‘dsyevr’ (default), ‘dsyev’, ‘dsyevd’, ‘dsyevx’
jobz [character(len=1)] – optional: ‘V’ (default) eigenvectors too, ‘N’ eigenvalues only
uplo [character(len=1)] – optional: ‘L’ (default) or ‘U’, triangle of A to be used
vl [real] – optional: lower bound of the eigenvalue interval (vl,vu]
vu [real] – optional: upper bound of the eigenvalue interval (vl,vu]
il [integer] – optional: index of the first eigenvalue wanted
iu [integer] – optional: index of the last eigenvalue wanted
tol [real] – optional: accepted but currently not used

interface sf_linalg/eigh_jacobi(a, d, v, nrot, u, sweep)

This subroutine computes all the eigenvalues and eigenvectors of a real symmetric (d_jacobi) or complex hermitian (c_jacobi) matrix with the Jacobi method of plane rotations. Only the upper triangle of the input matrix is referenced and the matrix is destroyed. The eigenvalues are NOT sorted.

real case: the eigenvectors are stored column-wise in v, v(:,j) being associated to d(j). The number of rotations performed is returned in nrot
complex case: the eigenvectors are stored in the ROWS of U, conjugated, so that \(A = U^\dagger \mathrm{diag}(D) U\). The sweep at which convergence was detected is returned in sweep

The program stops if the matrix is not square or if convergence is not reached within 50 sweeps.

Parameters:

a (•, •) [real/complex, inout] – in: symmetric matrix (upper triangle used); destroyed on exit
d (size(a, 1)) [real, out] – eigenvalues, NOT sorted, size n
v (size(a, 1), size(a, 2)) [real, out] – eigenvectors, v(:,j) <-> d(j), [n,n]
nrot [integer, out] – number of Jacobi rotations performed
u (size(a, 1), size(a, 2)) [complex, out] – eigenvectors: conjg(U(j,:)) <-> D(j), [n,n]
sweep [integer] – sweep at which convergence was detected

interface sf_linalg/eigvals(a)

This function returns the eigenvalues of a general real or complex square matrix A, computed with the LAPACK routines dgeev (deigvals) and zgeev (zeigvals) without computing the eigenvectors. The input matrix is not modified. The result is a complex array of size n, in the order given by LAPACK (not sorted). The program stops if A is not square or if the LAPACK routine returns an error code.

Parameters:: a (•, •) [real/complex, in] – eigenvalues of A, A c = lam c, size n
Result:: lam (size(a, 1)) [complex]

interface sf_linalg/eigvalsh(a)

This function returns the eigenvalues of a real symmetric (deigvalsh) or complex hermitian (zeigvalsh) matrix A, computed with the LAPACK routines dsyevd and zheevd without computing the eigenvectors. Only the upper triangle of A is referenced and the input matrix is not modified. The result is a real array, in ascending order. The program stops if A is not square or if the LAPACK routine returns an error code.

Parameters:: a (•, •) [real/complex, in] – eigenvalues of A, ascending, size n
Result:: lam (size(a, 1)) [real]

interface sf_linalg/svdvals(a)

This function returns the singular values of a real (dsvdvals) or complex (zsvdvals) \(m\times n\) matrix A, in descending order, computed with the LAPACK routine dgesvd (zgesvd) without computing the singular vectors. The input matrix is not modified. The result is a real array allocated with size min(m,n). The program stops if the LAPACK routine returns an error code.

Parameters:: a (•, •) [real/complex, in] – real matrix [m,n], not modified
Result:: s (•) [real, allocatable] – = lda

interface sf_linalg/svd(a, s, u, vtransp)

This subroutine computes the singular value decomposition \(A = U \Sigma V^T\) of a real (dsvd) or \(A = U \Sigma V^H\) of a complex (zsvd) \(m\times n\) matrix A, with the LAPACK routines dgesvd and zgesvd. The input matrix is not modified.

The singular values s are in descending order, and \(\Sigma\) is the \(m\times n\) matrix with \(\Sigma_{ii} = s_i\). The matrices U and Vtransp are full, with shapes [m,m] and [n,n]. Note that the routine returns the transpose (conjugate transpose in the complex case) of \(V\), not \(V\). The program stops if U or Vtransp do not have these shapes, or if the LAPACK routine returns an error code.

Parameters:

a (•, •) [real/complex, in] – real matrix [m,n], not modified
s (•) [real, out] – = lda
u (•, •) [real/complex, out] – left singular vectors (columns), [m,m]
vtransp (•, •) [real/complex, out] – V^T (V^H in the complex case), [n,n]

interface sf_linalg/inv(am)

This subroutine inverts in place a general real (dinv) or complex (zinv) square matrix Am using the LU factorization with partial pivoting, dgetrf + dgetri (zgetrf + zgetri). On output Am is replaced by its inverse. The program stops if Am is not square, if it is singular, or if the LAPACK routines return an error code.

Parameters:: am (•, •) [real/complex, inout] – in: square matrix to be inverted; out: its inverse

interface sf_linalg/p_inv(a, nblock)

This subroutine is the distributed-memory counterpart of inv, for real (p_dinv) and complex (p_zinv) square matrices, using ScaLAPACK. It is available only when compiling with _SCALAPACK.

The matrix A, a full copy of which is expected on every process, is distributed on the BLACS process grid of the module in a block-cyclic layout with square blocks of size Nblock, and inverted in place with the LU factorization, PDGETRF + PDGETRI (PZGETRF + PZGETRI). The inverse is gathered back into A. The program stops if A is not square or if a ScaLAPACK routine fails.

Parameters:

a (•, •) [real/complex, inout] – in: square matrix to be inverted; out: its inverse
nblock [integer] – block size of the block-cyclic distribution

interface sf_linalg/inv_sym(a[, uplo])

This subroutine inverts in place a real symmetric (dinv_sym) or complex symmetric (zinv_sym, NOT hermitian, see inv_her) matrix A using the Bunch-Kaufman factorization, dsytrf + dsytri (zsytrf + zsytri). Only the triangle selected by uplo (default 'U') is referenced on input. On output A contains the full inverse, the computed triangle being mirrored onto the other one. The program stops if a LAPACK routine returns an error code.

Parameters:: a (•, •) [real/complex] – in: symmetric matrix; out: its inverse (full matrix)
Options:: uplo [character(len=*)] – optional: ‘U’ (default) or ‘L’, triangle of A referenced

interface sf_linalg/inv_her(a[, uplo])

This subroutine inverts in place a complex hermitian matrix A using the Bunch-Kaufman factorization, zhetrf + zhetri. The hermiticity of A is tested first. Only the triangle selected by uplo (default 'U') is referenced by LAPACK, and on output A contains the full inverse, the other triangle being the complex conjugate of the computed one. The program stops if A is not hermitian or if a LAPACK routine returns an error code.

Parameters:: a (•, •) [complex] – in: hermitian matrix; out: its inverse (full matrix)
Options:: uplo [character(len=*)] – optional: ‘U’ (default) or ‘L’, triangle of A referenced

interface sf_linalg/inv_triang(a[, uplo, diag])

This subroutine inverts in place a real (dinv_triang) or complex (zinv_triang) triangular matrix A using the LAPACK routine dtrtri (ztrtri). On output the triangle selected by uplo (default 'U') contains the triangle of the inverse, the other one is not referenced. If diag='U' the matrix is assumed unit triangular and its diagonal is not referenced, the default being 'N'. The program stops if the LAPACK routine returns an error code, e.g. for a singular matrix.

Parameters:

a (•, •) [real/complex] – in: triangular matrix [n,n]; out: its inverse (same triangle)

Options:

uplo [character(len=1)] – optional: ‘U’ (default) upper, ‘L’ lower triangular
diag [character(len=1)] – optional: ‘N’ (default) non-unit, ‘U’ unit triangular

interface sf_linalg/inv_gj(a)

This subroutine inverts in place a real (dinv_gj) or complex (zinv_gj) square matrix a with the Gauss-Jordan elimination with full pivoting. On output a is replaced by its inverse. The program stops if the matrix is singular.

Parameters:: a (•, •) [real/complex, inout] – in: square matrix [n,n]; out: its inverse

interface sf_linalg/solve(a, b[, trans])

This subroutine solves the real or complex linear system \(A x = b\) using the LU factorization, dgetrf + dgetrs (zgetrf + zgetrs). The specific procedures cover a single right-hand side (dsolve_1rhs, zsolve_1rhs) and several right-hand sides stored as the columns of b (dsolve_mrhs, zsolve_mrhs).

The solution overwrites b. Note that A is factorized in place: despite the intent(in) it contains the LU factors on exit, so a copy must be passed if the original matrix is needed afterwards. If trans is 'T' or 'C' the transposed system \(A^T x = b\) is solved instead, the default being 'N'. The program stops if a LAPACK routine returns an error code, e.g. for a singular matrix.

Parameters:

a (•, •) [real/complex, in] – square matrix [n,n]; overwritten by its LU factors
b (various shapes) [real/complex, inout] – in: right-hand side, size n; out: solution x

Options:

trans [character(len=1)] – optional: ‘N’ (default) solve A x = b, ‘T’/’C’ solve A^T x = b

interface sf_linalg/lstsq(a, b)

This function returns the least-squares solution \(x = \mathrm{argmin}_x \| A x - b \|_2\) of the real (dlstsq) or complex (zlstsq) linear system \(A x = b\), with \(A\) an \(m\times n\) matrix of possibly non-square shape. It uses the LAPACK routine dgelsy (zgelsy), a complete orthogonal factorization with column pivoting, which is suitable also for rank-deficient matrices. The input arrays are not modified, and the solution is allocated with size n. The program stops if the LAPACK routine returns an error code.

Parameters:

a (•, •) [real/complex, in] – real matrix [m,n], not modified
b (•) [real/complex, in] – right-hand side, size m, not modified

Result:

x (•) [real/complex, allocatable] – least-squares solution, size n

interface sf_linalg/inv_tridiag(sub_, diag_, over_, inv, ainv, amat, nb, n[, n, nb])

This subroutine returns the diagonal elements, or the diagonal blocks, of the inverse of a real or complex (block) tridiagonal matrix, without building the full inverse, using a downward and an upward recursion. The specific procedures cover

a tridiagonal matrix given by its sub-diagonal sub_, main diagonal diag_ and over-diagonal over_, of size N: the N diagonal elements of the inverse are returned in Inv
a block tridiagonal matrix given by its block diagonals, made of Nb blocks of size \(N\times N\): the Nb diagonal blocks of the inverse are returned in Ainv
a full tridiagonal matrix Amat, or a full block tridiagonal matrix Amat made of Nb blocks of size \(N\times N\), treated in place: the diagonals are extracted with get_tridiag and on output Amat is a diagonal (block diagonal) matrix containing the diagonal elements (blocks) of the inverse

The tridiagonal structure of the input is not checked, see check_tridiag. The program stops if a zero pivot is found, i.e. if the matrix is ill-conditioned, or if Amat is not square.

Options:

n [integer, optional,default=1 + shape(sub, 0)] – size of the matrix
nb [integer, optional,default=1 + shape(sub, 0)] – number of blocks

Parameters:

sub_ (various shapes) [real/complex] – sub-diagonal, A(i+1,i), size N-1
diag_ (various shapes) [real/complex] – main diagonal, size N
over_ (various shapes) [real/complex] – over-diagonal, A(i,i+1), size N-1
inv (n) [real/complex] – diagonal elements of the inverse matrix, size N
ainv (nb, n, n) [real/complex] – diagonal blocks of the inverse matrix, [Nb,N,N]
amat (various shapes) [real/complex, inout] – in: tridiagonal [N,N]; out: diag. matrix of inverse diagonal

interface sf_linalg/check_tridiag(amat, nblock, nsize)

This function checks whether a real or complex matrix Amat is tridiagonal, i.e. whether all the elements outside the main, sub- and over-diagonals are exactly zero. The specific procedures cover a plain matrix and a matrix made of Nblock blocks of size \(N_{size}\times N_{size}\), to be checked for block tridiagonal structure. The result is .true. for a (block) tridiagonal matrix. The program stops if a plain matrix is not square.

Parameters:

amat (various shapes) [real/complex] – square matrix [N,N] to be checked
nblock [integer, in,required] – number of blocks
nsize [integer, in,required] – size of each block

Result:

mcheck [logical] – .true. if Amat is tridiagonal

interface sf_linalg/get_tridiag(amat, sub, diag, nblock, nsize[, over])

This subroutine extracts the three diagonals of a real or complex tridiagonal matrix Amat: the sub-diagonal sub(i) = A(i+1,i), the main diagonal diag(i) = A(i,i) and, optionally, the over-diagonal over(i) = A(i,i+1). The block version, for a matrix made of Nblock blocks of size \(N_{size}\times N_{size}\), extracts the block diagonals sub(i,:,:) = A(i+1,i), diag(i,:,:) = A(i,i) and over(i,:,:) = A(i,i+1). The input matrix is not modified and its structure is not checked, see check_tridiag. The inverse operation is build_tridiag.

Parameters:

amat (various shapes) [real/complex, in,required] – tridiagonal matrix [N,N]
sub (various shapes) [real/complex] – sub-diagonal, A(i+1,i), size N-1
diag (various shapes) [real/complex] – main diagonal, size N
nblock [integer, in,required,default=1 + shape(sub, 0)] – number of blocks
nsize [integer, in,required] – size of each block

Options:

over (various shapes) [real/complex] – optional out: over-diagonal, A(i,i+1), size N-1

interface sf_linalg/build_tridiag(sub, diag[, over, nblock, nsize])

This function builds a real or complex tridiagonal matrix from its three diagonals: the sub-diagonal A(i+1,i) = sub(i), the main diagonal A(i,i) = diag(i) and the over-diagonal A(i,i+1) = over(i). If over is not present the matrix is symmetric, over = sub. The block version, with Nblock blocks of size \(N_{size}\times N_{size}\), builds a block tridiagonal matrix from the block diagonals sub(i,:,:), diag(i,:,:) and over(i,:,:). The inverse operation is get_tridiag.

Parameters:

sub (various shapes) [real/complex] – sub-diagonal, size N-1
diag (various shapes) [real/complex, in,required] – main diagonal, size N

Options:

over (various shapes) [real/complex] – optional: over-diagonal, size N-1 (default: over = sub)
nblock [integer, optional,default=1 + shape(sub, 0)] – number of blocks
nsize [integer] – size of each block

Result:

amat (various shapes) [real/complex] – tridiagonal matrix, [N,N]

interface sf_linalg/det(a)

This function returns the determinant of a real (ddet) or complex (zdet) square matrix A, computed from its LU factorization (dgetrf, zgetrf) as the product of the diagonal elements of \(U\) with the sign changes due to the row permutations. The input matrix is not modified. The program stops if A is not square, or if the LAPACK routine returns an error code, which includes the case of an exactly singular matrix.

Parameters:: a (•, •) [real/complex, in] – square matrix [n,n], not modified
Result:: x [real/complex] – determinant of A

interface sf_linalg/deye(n, i, j)

This function returns real identity matrices or matrix elements: if one integer input variable is provided, it returns the the \(n\times n\) identity matrix deye(n) (deye_matrix). If instead two integer input variables are provided, it returns the Kronecker delta \(\delta_{ij}\), 1 if i==j and 0 otherwise (deye_indices). It is the same as eye, and the complex version is zeye.

Parameters:

n [integer, in] – size of the matrix
i [integer, in] – first index
j [integer, in] – second index

Result:

a (various shapes) [real] – identity matrix [n,n] or matrix element (scalar)

interface sf_linalg/zeye(n, i, j)

This function returns complex identity matrices or matrix elements: if one integer input variable is provided, it returns the the \(n\times n\) identity matrix zeye(n) (zeye_matrix). If instead two integer input variables are provided, it returns the Kronecker delta \(\delta_{ij}\), 1 if i==j and 0 otherwise (zeye_indices).

Parameters:

n [integer, in] – size of the matrix
i [integer, in] – first index
j [integer, in] – second index

Result:

a (various shapes) [complex] – identity matrix [n,n] or matrix element (scalar)

interface sf_linalg/eye(n, i, j)

This function returns real identity matrices or matrix elements: if one integer input variable is provided, it returns the the \(n\times n\) identity matrix deye(n) (deye_matrix). If instead two integer input variables are provided, it returns the Kronecker delta \(\delta_{ij}\), 1 if i==j and 0 otherwise (deye_indices). It is the same as deye, and the complex version is zeye.

Parameters:

n [integer, in] – size of the matrix
i [integer, in] – first index
j [integer, in] – second index

Result:

a (various shapes) [real] – identity matrix [n,n] or matrix element (scalar)

interface sf_linalg/diag(x)

This function builds a real (ddiag) or complex (zdiag) square matrix with the elements of the array x on the main diagonal and zeros elsewhere. The inverse operation is diagonal.

Parameters:: x (•) [real/complex, in] – diagonal elements, size n
Result:: a (•, •) [real/complex, allocatable] – diagonal matrix, A(i,i) = x(i), [n,n]

interface sf_linalg/diagonal(a)

This function returns the main diagonal of a real (d_diagonal) or complex (z_diagonal) matrix A, expected to be square: the size of the result is size(A,1). The inverse operation is diag.

Parameters:: a (•, •) [real/complex, in] – square matrix [n,n]
Result:: dd (size(a, 1)) [real/complex] – diagonal elements, dd(i) = A(i,i), size n

interface sf_linalg/trace(a)

This function returns the trace of a real (dtrace) or complex (ztrace) matrix A, i.e. the sum of the elements along the main diagonal. For a non-square matrix the sum runs up to min(m,n).

Parameters:: a (•, •) [real/complex, in] – real matrix [m,n]
Result:: t [real/complex] – trace of A

interface sf_linalg/zeros(n, n1, n2, n3, n4, n5, n6, n7)

This function returns a complex(8) array with all the elements set to 0. The shape is given by one integer extent for each dimension, zeros(n1,...,nk), and the ranks from 1 to 7 are covered.

Parameters:

n [integer, in] – extent of dimension 1
n1 [integer, in] – extent of dimension 1
n2 [integer, in] – extent of dimension 2
n3 [integer, in] – extent of dimension 3
n4 [integer, in] – extent of dimension 4
n5 [integer, in] – extent of dimension 5
n6 [integer, in] – extent of dimension 6
n7 [integer, in] – extent of dimension 7

Result:

a (various shapes) [complex] – complex array, [n]

interface sf_linalg/ones(n, n1, n2, n3, n4, n5, n6, n7)

This function returns a complex(8) array with all the elements set to 1. The shape is given by one integer extent for each dimension, ones(n1,...,nk), and the ranks from 1 to 7 are covered.

Parameters:

n [integer, in] – extent of dimension 1
n1 [integer, in] – extent of dimension 1
n2 [integer, in] – extent of dimension 2
n3 [integer, in] – extent of dimension 3
n4 [integer, in] – extent of dimension 4
n5 [integer, in] – extent of dimension 5
n6 [integer, in] – extent of dimension 6
n7 [integer, in] – extent of dimension 7

Result:

a (various shapes) [complex] – complex array, [n]

interface sf_linalg/kron(a, b)

This function returns the Kronecker (tensor) product \(A \otimes B\) of two matrices A and B, i.e. the block matrix whose \((i,j)\) block is \(A_{ij} B\). The result has shape [size(A,1)*size(B,1),size(A,2)*size(B,2)]. The specific procedures cover integer, real, real \(\otimes\) complex, complex \(\otimes\) real and complex matrices, the result being complex if any of the two is complex. It is the same as kronecker_product and as the operator .kx..

Parameters:

a (•, •) [integer/real/complex, in] – first factor, [rowA,colA]
b (•, •) [integer/real/complex, in] – Kronecker product, [rowA*rowB,colA*colB]

Result:

axb (size(a, 1) · size(b, 1), size(a, 2) · size(b, 2)) [integer/real/complex]

interface sf_linalg/kronecker_product(a, b)

This function returns the Kronecker (tensor) product \(A \otimes B\) of two matrices A and B, i.e. the block matrix whose \((i,j)\) block is \(A_{ij} B\). The result has shape [size(A,1)*size(B,1),size(A,2)*size(B,2)]. The specific procedures cover integer, real, real \(\otimes\) complex, complex \(\otimes\) real and complex matrices, the result being complex if any of the two is complex. It is the same as kron and as the operator .kx..

Parameters:

a (•, •) [integer/real/complex, in] – first factor, [rowA,colA]
b (•, •) [integer/real/complex, in] – Kronecker product, [rowA*rowB,colA*colB]

Result:

axb (size(a, 1) · size(b, 1), size(a, 2) · size(b, 2)) [integer/real/complex]

interface sf_linalg/operator(.kx.)(a, b)

This operator returns the Kronecker (tensor) product \(A \otimes B\) of two matrices A and B, i.e. the block matrix whose \((i,j)\) block is \(A_{ij} B\). The result has shape [size(A,1)*size(B,1),size(A,2)*size(B,2)]. The specific procedures cover integer, real, real \(\otimes\) complex, complex \(\otimes\) real and complex matrices, the result being complex if any of the two is complex. It is used as AxB = A .kx. B, and it is the same as kron.

Parameters:

a (•, •) [integer/real/complex, in] – first factor, [rowA,colA]
b (•, •) [integer/real/complex, in] – Kronecker product, [rowA*rowB,colA*colB]

Result:

axb (size(a, 1) · size(b, 1), size(a, 2) · size(b, 2)) [integer/real/complex]

interface sf_linalg/outerprod(a, b)

This function returns the outer product of two real or complex one-dimensional arrays a and b, i.e. the matrix \(A_{ij} = a_i b_j\), of shape [size(a),size(b)]. No complex conjugation is applied.

Parameters:

a (•) [real/complex, in]
b (•) [real/complex, in]

Result:

outerprod (size(a), size(b)) [real/complex]

interface sf_linalg/cross_product(a, b)

This function returns the cross (vector) product \(c = a \times b\) of two real or complex vectors a and b with 3 components, \(c_1 = a_2 b_3 - a_3 b_2\), \(c_2 = a_3 b_1 - a_1 b_3\), \(c_3 = a_1 b_2 - a_2 b_1\). Only the 3-dimensional version is available through this interface.

Parameters:

a (3) [real/complex]
b (3) [real/complex]

Result:

c (3) [real/complex]

interface sf_linalg/s3_product(a, b, c)

This function returns the triple product \(a \cdot (b \times c)\) of three real or complex vectors a, b and c with 3 components, evaluated as dot_product(a,cross_product(b,c)). The result is a real number: in the complex case a is conjugated by dot_product and the imaginary part of the result is discarded.

Parameters:

a (3) [real/complex, in]
b (3) [real/complex, in]
c (3) [real/complex, in]

Result:

s3 [real]

interface sf_linalg/mat_product(a, b, c[, alfa, beta])

This subroutine computes the real (d_matmul) or complex (z_matmul) matrix-matrix product \(C = \alpha A B + \beta C\) with the BLAS routine dgemm (zgemm), \(A\) being \(N\times K\), \(B\) \(K\times M\) and \(C\) \(N\times M\). By default \(\alpha=1\) and \(\beta=0\), so that C contains the plain product \(A B\) on output. A and B are not modified. The program stops if B or C do not have the shapes [K,M] and [N,M]. The function form is the operator .x..

Parameters:

a (•, •) [real/complex, inout] – matrix [N,K], not modified
b (•, •) [real/complex, inout] – matrix [K,M], not modified
c (•, •) [real/complex, inout] – in/out: C = alfa*A*B + beta*C, [N,M]

Options:

alfa [real/complex] – optional scalar, default 1
beta [real/complex] – optional scalar, default 0

interface sf_linalg/operator(.x.)(a, b)

This operator returns the real or complex matrix-matrix product \(C = A B\), computed with the BLAS routine dgemm (zgemm), \(A\) being \(N\times K\) and \(B\) \(K\times M\). It is used as C = A .x. B, through the functions d_matmul_ and z_matmul_. The program stops if B does not have the shape [K,M]. The subroutine form is mat_product.

Parameters:

a (•, •) [real/complex, in] – matrix [N,K]
b (•, •) [real/complex, in] – matrix [K,M]

Result:

c (size(a, 1), size(b, 2)) [real/complex] – product A*B, [N,M]

interface sf_linalg/p_mat_product(a, b, c, nblock[, alfa, beta])

This subroutine is the distributed-memory counterpart of mat_product: it computes the real (p_d_matmul) or complex (p_z_matmul) matrix-matrix product \(C = \alpha A B + \beta C\) with the ScaLAPACK routine PDGEMM (PZGEMM). It is available only when compiling with _SCALAPACK.

The matrices A and B, full copies of which are expected on every process, are distributed on the BLACS process grid of the module in a block-cyclic layout with square blocks of size Nblock, and the product is gathered back into C. By default \(\alpha=1\) and \(\beta=0\). The program stops if B or C do not have the shapes [K,M] and [N,M]. The function form is the operator .Px..

Parameters:

a (•, •) [real/complex, inout] – matrix [N,K]
b (•, •) [real/complex, inout] – matrix [K,M]
c (various shapes) [real/complex, inout] – out: product alfa*A*B, [N,M]
nblock [integer] – block size of the block-cyclic distribution

Options:

alfa [real] – optional scalar, default 1
beta [real] – optional scalar, default 0

interface sf_linalg/operator(.px.)(a, b)

This operator returns the real or complex matrix-matrix product \(C = A B\), computed in parallel with the ScaLAPACK routine PDGEMM (PZGEMM), \(A\) being \(N\times K\) and \(B\) \(K\times M\). It is used as C = A .Px. B, through the functions p_d_matmul_f and p_z_matmul_f, and it is available only when compiling with _SCALAPACK. The block size of the block-cyclic distribution is chosen automatically, as the largest power of 2 smaller than min(N,K,M), up to 64. The program stops if B does not have the shape [K,M]. The subroutine form is p_mat_product.

Parameters:

a (•, •) [real/complex, in] – matrix [N,K]
b (•, •) [real/complex, in] – matrix [K,M]

Result:

c (size(a, 1), size(b, 2)) [real/complex] – product A*B, [N,M]

function sf_linalg/deye_tridiag(nblock, n)

This function returns Nblock real \(N\times N\) identity matrices stacked in a rank-3 array, i.e. the blocks of a block identity matrix. It can be used to build block tridiagonal matrices, see build_tridiag.

Parameters:

nblock [integer] – number of blocks
n [integer] – size of each block

Result:

eye_block (nblock, n, n) [real] – identity blocks, eye_block(i,:,:) = 1, [Nblock,N,N]

function sf_linalg/zeye_tridiag(nblock, n)

This function is the complex counterpart of deye_tridiag: it returns Nblock complex \(N\times N\) identity matrices stacked in a rank-3 array.

Parameters:

nblock [integer] – number of blocks
n [integer] – size of each block

Result:

eye_block (nblock, n, n) [complex] – identity blocks, eye_block(i,:,:) = 1, [Nblock,N,N]