The stdlib linear algebra library provides high-level APIs for dealing with common linear algebra operations.
Experimental
BLAS and LAPACK backends provide efficient low level implementations of many linear algebra algorithms, and are employed for non-trivial operators.
A Modern Fortran version of the Reference-LAPACK 3.10.1 implementation is provided as a backend.
Modern Fortran modules with full explicit typing features are provided after an
automated conversion
of the legacy codes:
- [stdlib_linalg_blas(module)], [stdlib_linalg_lapack(module)] provide kind-agnostic interfaces to all functions.
- Both libraries are available for 32- (sp), 64- (dp) and 128-bit (qp) real and complex numbers (the latter if available in the current build)
- Free format, lower-case style
- implicit none(type, external) applied to all procedures and modules
- intent added and all pure procedures where possible
- stdlib provides all procedures in two different flavors: (a) original BLAS/LAPACK names with a prefix stdlib_?<name> (ex: stdlib_dgemv, stdlib_sgemv); (b) A generic, kind agnostic <name>, i.e. gemv.
- F77-style parameters removed, and all numeric constants have been generalized with KIND-dependent Fortran intrinsics.
- preprocessor-based OpenMP directives retained.
The single-source module structure hopefully allows for cross-procedural inlining which is otherwise impossible without link-time optimization.
When available, highly optimized libraries that take advantage of specialized processor instructions should be preferred over the stdlib implementation.
Examples of such libraries are: OpenBLAS, MKL (TM), Accelerate, and ATLAS. In order to enable their usage, simply ensure that the following pre-processor macros are defined:
STDLIB_EXTERNAL_BLAS wraps all BLAS procedures (except for the 128-bit ones) to an external librarySTDLIB_EXTERNAL_LAPACK wraps all LAPACK procedures (except for the 128-bit ones) to an external libraryThese can be enabled during the build process. For example, with CMake, one can enable these preprocessor directives using add_compile_definitions(STDLIB_EXTERNAL_BLAS STDLIB_EXTERNAL_LAPACK).
The same is possible from the fpm branch, where the cpp preprocessor is enabled by default. For example, the macros can be added to the project's manifest:
# Link against appropriate external BLAS and LAPACK libraries, if necessary
[build]
link = ["blas", "lapack"]
[dependencies]
stdlib="*"
# Macros are only needed if using an external library
[preprocess]
[preprocess.cpp]
macros = ["STDLIB_EXTERNAL_BLAS", "STDLIB_EXTERNAL_LAPACK"]
or directly via compiler flags:
fpm build --flag "-DSTDLIB_EXTERNAL_BLAS -DSTDLIB_EXTERNAL_LAPACK -lblas -llapack".
All procedures in the BLAS and LAPACK backends follow the standard interfaces from the
Reference LAPACK. So, the online Users Guide
should be consulted for the full API and descriptions of procedure arguments and their usage.
The stdlib implementation makes both kind-agnostic and specific procedure interfaces available via modules
[stdlib_linalg_blas(module)] and [stdlib_linalg_lapack(module)]. Because all procedures start with a letter
that indicates the base datatype, the stdlib generic
interface drops the heading letter and contains all kind-dependent implementations. For example, the generic
interface to the axpy function looks like:
!> AXPY: constant times a vector plus a vector.
interface axpy
module procedure stdlib_saxpy
module procedure stdlib_daxpy
module procedure stdlib_qaxpy
module procedure stdlib_caxpy
module procedure stdlib_zaxpy
module procedure stdlib_waxpy
end interface axpy
The generic interface is the endpoint for using an external library. Whenever the latter is used, references
to the internal module procedures are replaced with interfaces to the external library,
for example:
!> AXPY: constant times a vector plus a vector.
interface axpy
pure subroutine caxpy(n,ca,cx,incx,cy,incy)
import sp,dp,qp,ilp,lk
implicit none(type,external)
complex(sp), intent(in) :: ca,cx(*)
integer(ilp), intent(in) :: incx,incy,n
complex(sp), intent(inout) :: cy(*)
end subroutine caxpy
! [....]
module procedure stdlib_qaxpy
end interface axpy
Note that the 128-bit functions are only provided by stdlib and always point to the internal implementation.
Because 128-bit precision is identified as [stdlib_kinds(module):qp], initials for 128-bit procedures were
labelled as q (quadruple-precision reals) and w ("wide" or quadruple-precision complex numbers).
Extended precision ([stdlib_kinds(module):xdp]) calculations are currently not supported.
program example_gemv
use stdlib_linalg, only: eye
use stdlib_linalg_blas, only: sp,gemv
implicit none(type,external)
real(sp) :: A(2, 2), B(2), C(2)
B = [1.0,2.0]
A = eye(2)
! Use legacy BLAS interface
call gemv('No transpose',m=size(A,1),n=size(A,2),alpha=1.0,a=A,lda=size(A,1),x=B,incx=1,beta=0.0,y=C,incy=1)
print *, C ! returns 1.0 2.0
end program example_gemv
program example_getrf
use stdlib_linalg, only: eye
use stdlib_linalg_lapack, only: dp,ilp,getrf
implicit none(type,external)
real(dp) :: A(3, 3)
integer(ilp) :: ipiv(3),info
A = eye(3)
! LAPACK matrix factorization interface (overwrite result)
call getrf(size(A,1),size(A,2),A,size(A,1),ipiv,info)
print *, info ! info==0: Success!
end program example_getrf
The Fortran Standard Library is distributed under the MIT License. LAPACK and its contained BLAS are a
freely-available software package. They are available from netlib via anonymous
ftp and the World Wide Web. Thus, they can be included in commercial software packages (and have been).
The license used for the BLAS and LAPACK backends is the modified BSD license.
The header of the LICENSE.txt file has as its licensing requirements:
Copyright (c) 1992-2013 The University of Tennessee and The University
of Tennessee Research Foundation. All rights
reserved.
Copyright (c) 2000-2013 The University of California Berkeley. All
rights reserved.
Copyright (c) 2006-2013 The University of Colorado Denver. All rights
reserved.
$COPYRIGHT$
Additional copyrights may follow
$HEADER$
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
- Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer listed
in this license in the documentation and/or other materials
provided with the distribution.
- Neither the name of the copyright holders nor the names of its
contributors may be used to endorse or promote products derived from
this software without specific prior written permission.
The copyright holders provide no reassurances that the source code
provided does not infringe any patent, copyright, or any other
intellectual property rights of third parties. The copyright holders
disclaim any liability to any recipient for claims brought against
recipient by any third party for infringement of that parties
intellectual property rights.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
So the license for the LICENSE.txt code is compatible with the use of
modified versions of the code in the Fortran Standard Library under the MIT license.
Credit for the BLAS, LAPACK libraries should be given to the LAPACK authors.
According to the original license, we also changed the name of the routines and commented the changes made
to the original.
diag - Create a diagonal array or extract the diagonal elements of an arrayExperimental
Create a diagonal array or extract the diagonal elements of an array
d = diag (a [, k])
a: Shall be a rank-1 or or rank-2 array. If a is a rank-1 array (i.e. a vector) then diag returns a rank-2 array with the elements of a on the diagonal. If a is a rank-2 array (i.e. a matrix) then diag returns a rank-1 array of the diagonal elements.
k (optional): Shall be a scalar of type integer and specifies the diagonal. The default k = 0 represents the main diagonal, k > 0 are diagonals above the main diagonal, k < 0 are diagonals below the main diagonal.
Returns a diagonal array or a vector with the extracted diagonal elements.
program example_diag1
use stdlib_linalg, only: diag
implicit none
real, allocatable :: A(:, :)
integer :: i
A = diag([(1, i=1, 10)]) ! creates a 10 by 10 identity matrix
end program example_diag1
program example_diag2
use stdlib_linalg, only: diag
implicit none
real, allocatable :: v(:)
real, allocatable :: A(:, :)
v = [1, 2, 3, 4, 5]
A = diag(v) ! creates a 5 by 5 matrix with elements of v on the diagonal
end program example_diag2
program example_diag3
use stdlib_linalg, only: diag
implicit none
integer, parameter :: n = 10
real :: c(n), ul(n - 1)
real :: A(n, n)
c = 2
ul = -1
A = diag(ul, -1) + diag(c) + diag(ul, 1) ! Gil Strang's favorite matrix
end program example_diag3
program example_diag4
use stdlib_linalg, only: diag
implicit none
integer, parameter :: n = 12
real :: A(n, n)
real :: v(n)
call random_number(A)
v = diag(A) ! v contains diagonal elements of A
end program example_diag4
program example_diag5
use stdlib_linalg, only: diag
implicit none
integer, parameter :: n = 3
real :: A(n, n)
real, allocatable :: v(:)
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [n, n])
v = diag(A, -1) ! v is [2,6]
v = diag(A, 1) ! v is [4,8]
end program example_diag5
eye - Construct the identity matrixExperimental
Pure function.
Construct the identity matrix.
I = eye (dim1 [, dim2])
dim1: Shall be a scalar of default type integer.
This is an intent(in) argument.
dim2: Shall be a scalar of default type integer.
This is an intent(in) and optional argument.
Return the identity matrix, i.e. a matrix with ones on the main diagonal and zeros elsewhere. The return value is of type integer(int8).
The use of int8 was suggested to save storage.
Since the result of eye is of integer(int8) type, one should be careful about using it in arithmetic expressions. For example:
!> Be careful
A = eye(2,2)/2 !! A == 0.0
!> Recommend
A = eye(2,2)/2.0 !! A == diag([0.5, 0.5])
program example_eye1
use stdlib_linalg, only: eye
implicit none
integer :: i(2, 2)
real :: a(3, 3)
real :: b(2, 3) !! Matrix is non-square.
complex :: c(2, 2)
I = eye(2) !! [1,0; 0,1]
A = eye(3) !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
A = eye(3, 3) !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
B = eye(2, 3) !! [1.0,0.0,0.0; 0.0,1.0,0.0]
C = eye(2, 2) !! [(1.0,0.0),(0.0,0.0); (0.0,0.0),(1.0,0.0)]
C = (1.0, 1.0)*eye(2, 2) !! [(1.0,1.0),(0.0,0.0); (0.0,0.0),(1.0,1.0)]
end program example_eye1
program example_eye2
use stdlib_linalg, only: eye, diag
implicit none
print *, all(eye(4) == diag([1, 1, 1, 1])) ! prints .true.
end program example_eye2
trace - Trace of a matrixExperimental
Trace of a matrix (rank-2 array)
result = trace (A)
A: Shall be a rank-2 array. If A is not square, then trace(A) will return the sum of diagonal values from the square sub-section of A.
Returns the trace of the matrix, i.e. the sum of diagonal elements.
program example_trace
use stdlib_linalg, only: trace
implicit none
real :: A(3, 3)
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
print *, trace(A) ! 1 + 5 + 9
end program example_trace
outer_product - Computes the outer product of two vectorsExperimental
Computes the outer product of two vectors
d = outer_product (u, v)
u: Shall be a rank-1 array
v: Shall be a rank-1 array
Returns a rank-2 array equal to u v^T (where u, v are considered column vectors). The shape of the returned array is [size(u), size(v)].
program example_outer_product
use stdlib_linalg, only: outer_product
implicit none
real, allocatable :: A(:, :), u(:), v(:)
u = [1., 2., 3.]
v = [3., 4.]
A = outer_product(u, v)
!A = reshape([3., 6., 9., 4., 8., 12.], [3,2])
end program example_outer_product
kronecker_product - Computes the Kronecker product of two rank-2 arraysExperimental
Computes the Kronecker product of two rank-2 arrays
C = kronecker_product (A, B)
A: Shall be a rank-2 array with dimensions M1, N1
B: Shall be a rank-2 array with dimensions M2, N2
Returns a rank-2 array equal to A \otimes B. The shape of the returned array is [M1*M2, N1*N2].
program example_kronecker_product
use stdlib_linalg, only: kronecker_product
implicit none
integer, parameter :: m1 = 1, n1 = 2, m2 = 2, n2 = 3
integer :: i, j
real :: A(m1, n1), B(m2,n2)
real, allocatable :: C(:,:)
do j = 1, n1
do i = 1, m1
A(i,j) = i*j ! A = [1, 2]
end do
end do
do j = 1, n2
do i = 1, m2 ! B = [ 1, 2, 3 ]
B(i,j) = i*j ! [ 2, 4, 6 ]
end do
end do
C = kronecker_product(A, B)
! C = [ a(1,1) * B(:,:) | a(1,2) * B(:,:) ]
! or in other words,
! C = [ 1.00 2.00 3.00 2.00 4.00 6.00 ]
! [ 2.00 4.00 6.00 4.00 8.00 12.00 ]
end program example_kronecker_product
cross_product - Computes the cross product of two vectorsExperimental
Computes the cross product of two vectors
c = cross_product (a, b)
a: Shall be a rank-1 and size-3 array
b: Shall be a rank-1 and size-3 array
Returns a rank-1 and size-3 array which is perpendicular to both a and b.
program demo_cross_product
use stdlib_linalg, only: cross_product
implicit none
real :: a(3), b(3), c(3)
a = [1., 0., 0.]
b = [0., 1., 0.]
c = cross_product(a, b)
!c = [0., 0., 1.]
end program demo_cross_product
is_square - Checks if a matrix is squareExperimental
Checks if a matrix is square
d = is_square (A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is square, and .false. otherwise.
program example_is_square
use stdlib_linalg, only: is_square
implicit none
real :: A(2, 2), B(3, 2)
logical :: res
A = reshape([1., 2., 3., 4.], shape(A))
B = reshape([1., 2., 3., 4., 5., 6.], shape(B))
res = is_square(A) ! returns .true.
res = is_square(B) ! returns .false.
end program example_is_square
is_diagonal - Checks if a matrix is diagonalExperimental
Checks if a matrix is diagonal
d = is_diagonal (A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is diagonal, and .false. otherwise.
Note that nonsquare matrices may be diagonal, so long as a_ij = 0 when i /= j.
program example_is_diagonal
use stdlib_linalg, only: is_diagonal
implicit none
real :: A(2, 2), B(2, 2)
logical :: res
A = reshape([1., 0., 0., 4.], shape(A))
B = reshape([1., 0., 3., 4.], shape(B))
res = is_diagonal(A) ! returns .true.
res = is_diagonal(B) ! returns .false.
end program example_is_diagonal
is_symmetric - Checks if a matrix is symmetricExperimental
Checks if a matrix is symmetric
d = is_symmetric (A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is symmetric, and .false. otherwise.
program example_is_symmetric
use stdlib_linalg, only: is_symmetric
implicit none
real :: A(2, 2), B(2, 2)
logical :: res
A = reshape([1., 3., 3., 4.], shape(A))
B = reshape([1., 0., 3., 4.], shape(B))
res = is_symmetric(A) ! returns .true.
res = is_symmetric(B) ! returns .false.
end program example_is_symmetric
is_skew_symmetric - Checks if a matrix is skew-symmetricExperimental
Checks if a matrix is skew-symmetric
d = is_skew_symmetric (A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is skew-symmetric, and .false. otherwise.
program example_is_skew_symmetric
use stdlib_linalg, only: is_skew_symmetric
implicit none
real :: A(2, 2), B(2, 2)
logical :: res
A = reshape([0., -3., 3., 0.], shape(A))
B = reshape([0., 3., 3., 0.], shape(B))
res = is_skew_symmetric(A) ! returns .true.
res = is_skew_symmetric(B) ! returns .false.
end program example_is_skew_symmetric
is_hermitian - Checks if a matrix is HermitianExperimental
Checks if a matrix is Hermitian
d = is_hermitian (A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is Hermitian, and .false. otherwise.
program example_is_hermitian
use stdlib_linalg, only: is_hermitian
implicit none
complex :: A(2, 2), B(2, 2)
logical :: res
A = reshape([cmplx(1., 0.), cmplx(3., -1.), cmplx(3., 1.), cmplx(4., 0.)], shape(A))
B = reshape([cmplx(1., 0.), cmplx(3., 1.), cmplx(3., 1.), cmplx(4., 0.)], shape(B))
res = is_hermitian(A) ! returns .true.
res = is_hermitian(B) ! returns .false.
end program example_is_hermitian
is_triangular - Checks if a matrix is triangularExperimental
Checks if a matrix is triangular
d = is_triangular (A,uplo)
A: Shall be a rank-2 array
uplo: Shall be a single character from {'u','U','l','L'}
Returns a logical scalar that is .true. if the input matrix is the type of triangular specified by uplo (upper or lower), and .false. otherwise.
Note that the definition of triangular used in this implementation allows nonsquare matrices to be triangular.
Specifically, upper triangular matrices satisfy a_ij = 0 when j < i, and lower triangular matrices satisfy a_ij = 0 when j > i.
program example_is_triangular
use stdlib_linalg, only: is_triangular
implicit none
real :: A(3, 3), B(3, 3)
logical :: res
A = reshape([1., 0., 0., 4., 5., 0., 7., 8., 9.], shape(A))
B = reshape([1., 0., 3., 4., 5., 0., 7., 8., 9.], shape(B))
res = is_triangular(A, 'u') ! returns .true.
res = is_triangular(B, 'u') ! returns .false.
end program example_is_triangular
is_hessenberg - Checks if a matrix is hessenbergExperimental
Checks if a matrix is Hessenberg
d = is_hessenberg (A,uplo)
A: Shall be a rank-2 array
uplo: Shall be a single character from {'u','U','l','L'}
Returns a logical scalar that is .true. if the input matrix is the type of Hessenberg specified by uplo (upper or lower), and .false. otherwise.
Note that the definition of Hessenberg used in this implementation allows nonsquare matrices to be Hessenberg.
Specifically, upper Hessenberg matrices satisfy a_ij = 0 when j < i-1, and lower Hessenberg matrices satisfy a_ij = 0 when j > i+1.
program example_is_hessenberg
use stdlib_linalg, only: is_hessenberg
implicit none
real :: A(3, 3), B(3, 3)
logical :: res
A = reshape([1., 2., 0., 4., 5., 6., 7., 8., 9.], shape(A))
B = reshape([1., 2., 3., 4., 5., 6., 7., 8., 9.], shape(B))
res = is_hessenberg(A, 'u') ! returns .true.
res = is_hessenberg(B, 'u') ! returns .false.
end program example_is_hessenberg
det - Computes the determinant of a square matrixExperimental
This function computes the determinant of a real or complex square matrix.
This interface comes with a pure version det(a), and a non-pure version det(a,overwrite_a,err) that
allows for more expert control.
c = det (a [, overwrite_a, err])
a: Shall be a rank-2 square array
overwrite_a (optional): Shall be an input logical flag. if .true., input matrix a will be used as temporary storage and overwritten. This avoids internal data allocation.
This is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.
Returns a real scalar value of the same kind of a that represents the determinant of the matrix.
Raises LINALG_ERROR if the matrix is singular.
Raises LINALG_VALUE_ERROR if the matrix is non-square.
Exceptions are returned to the err argument if provided; an error stop is triggered otherwise.
program example_determinant
use stdlib_kinds, only: dp
use stdlib_linalg, only: det, linalg_state_type
implicit none
type(linalg_state_type) :: err
real(dp) :: d
! Compute determinate of a real matrix
d = det(reshape([real(dp)::1,2,3,4],[2,2]))
print *, d ! a*d-b*c = -2.0
end program example_determinant
.det. - Determinant operator of a square matrixExperimental
This operator returns the determinant of a real square matrix.
This interface is equivalent to the pure version of determinant det.
c = [[stdlib_linalg(module):operator(.det.)(interface)]] (a)
a: Shall be a rank-2 square array of any real or complex kinds. It is an intent(in) argument.
Returns a real scalar value that represents the determinnt of the matrix.
Raises LINALG_ERROR if the matrix is singular.
Raises LINALG_VALUE_ERROR if the matrix is non-square.
Exceptions trigger an error stop.
program example_determinant2
use stdlib_kinds, only: dp
use stdlib_linalg, only: operator(.det.)
implicit none
real(dp) :: d
! Compute determinate of a real matrix
d = .det.reshape([real(dp)::1,2,3,4],[2,2])
print *, d ! a*d-b*c = -2.0
end program example_determinant2