linalg

Linear Algebra

The stdlib linear algebra library provides high-level APIs for dealing with common linear algebra operations.

BLAS and LAPACK

Status

Experimental

Description

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 library
  • STDLIB_EXTERNAL_LAPACK wraps all LAPACK procedures (except for the 128-bit ones) to an external library

These 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".

Syntax

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.

Example

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

Licensing

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 array

Status

Experimental

Description

Create a diagonal array or extract the diagonal elements of an array

Syntax

d = diag (a [, k])

Arguments

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.

Return value

Returns a diagonal array or a vector with the extracted diagonal elements.

Example

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 matrix

Status

Experimental

Class

Pure function.

Description

Construct the identity matrix.

Syntax

I = eye (dim1 [, dim2])

Arguments

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 value

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.

Warning

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])

Example

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 matrix

Status

Experimental

Description

Trace of a matrix (rank-2 array)

Syntax

result = trace (A)

Arguments

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.

Return value

Returns the trace of the matrix, i.e. the sum of diagonal elements.

Example

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 vectors

Status

Experimental

Description

Computes the outer product of two vectors

Syntax

d = outer_product (u, v)

Arguments

u: Shall be a rank-1 array

v: Shall be a rank-1 array

Return value

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)].

Example

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 arrays

Status

Experimental

Description

Computes the Kronecker product of two rank-2 arrays

Syntax

C = kronecker_product (A, B)

Arguments

A: Shall be a rank-2 array with dimensions M1, N1

B: Shall be a rank-2 array with dimensions M2, N2

Return value

Returns a rank-2 array equal to A \otimes B. The shape of the returned array is [M1*M2, N1*N2].

Example

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 vectors

Status

Experimental

Description

Computes the cross product of two vectors

Syntax

c = cross_product (a, b)

Arguments

a: Shall be a rank-1 and size-3 array

b: Shall be a rank-1 and size-3 array

Return value

Returns a rank-1 and size-3 array which is perpendicular to both a and b.

Example

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 square

Status

Experimental

Description

Checks if a matrix is square

Syntax

d = is_square (A)

Arguments

A: Shall be a rank-2 array

Return value

Returns a logical scalar that is .true. if the input matrix is square, and .false. otherwise.

Example

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 diagonal

Status

Experimental

Description

Checks if a matrix is diagonal

Syntax

d = is_diagonal (A)

Arguments

A: Shall be a rank-2 array

Return value

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.

Example

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 symmetric

Status

Experimental

Description

Checks if a matrix is symmetric

Syntax

d = is_symmetric (A)

Arguments

A: Shall be a rank-2 array

Return value

Returns a logical scalar that is .true. if the input matrix is symmetric, and .false. otherwise.

Example

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-symmetric

Status

Experimental

Description

Checks if a matrix is skew-symmetric

Syntax

d = is_skew_symmetric (A)

Arguments

A: Shall be a rank-2 array

Return value

Returns a logical scalar that is .true. if the input matrix is skew-symmetric, and .false. otherwise.

Example

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 Hermitian

Status

Experimental

Description

Checks if a matrix is Hermitian

Syntax

d = is_hermitian (A)

Arguments

A: Shall be a rank-2 array

Return value

Returns a logical scalar that is .true. if the input matrix is Hermitian, and .false. otherwise.

Example

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 triangular

Status

Experimental

Description

Checks if a matrix is triangular

Syntax

d = is_triangular (A,uplo)

Arguments

A: Shall be a rank-2 array

uplo: Shall be a single character from {'u','U','l','L'}

Return value

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.

Example

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 hessenberg

Status

Experimental

Description

Checks if a matrix is Hessenberg

Syntax

d = is_hessenberg (A,uplo)

Arguments

A: Shall be a rank-2 array

uplo: Shall be a single character from {'u','U','l','L'}

Return value

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.

Example

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