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 labelled as x (extended-precision reals). and y (extended-precision complex numbers).

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

Stable

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

Stable

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

Stable

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

solve - Solves a linear matrix equation or a linear system of equations.

Status

Stable

Description

This function computes the solution to a linear matrix equation $A \cdot x = b$, where $A$ is a square, full-rank, real or complex matrix.

Result vector or array x returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. An error is returned if the matrix is rank-deficient or singular to working precision. The solver is based on LAPACK's *GESV backends.

Syntax

Pure interface:

x = solve (a, b)

Expert interface:

x = solve (a, b [, overwrite_a], err)

Arguments

a: Shall be a rank-2 real or complex square array containing the coefficient matrix. It is normally an intent(in) argument. If overwrite_a=.true., it is an intent(inout) argument and is destroyed by the call.

b: Shall be a rank-1 or rank-2 array of the same kind as a, containing the right-hand-side vector(s). It is an intent(in) argument.

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. The function is not pure if this argument is provided.

Return value

For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.

Raises LINALG_ERROR if the matrix is singular to working precision. Raises LINALG_VALUE_ERROR if the matrix and rhs vectors have invalid/incompatible sizes. If err is not present, exceptions trigger an error stop.

Example

program example_solve1
  use stdlib_linalg_constants, only: sp
  use stdlib_linalg, only: solve, linalg_state_type 
  implicit none

  real(sp), allocatable :: A(:,:),b(:),x(:)

  ! Solve a system of 3 linear equations: 
  !  4x + 3y + 2z =  25
  ! -2x + 2y + 3z = -10
  !  3x - 5y + 2z =  -4

  ! Note: Fortran is column-major! -> transpose 
  A = transpose(reshape([ 4, 3, 2, &
                         -2, 2, 3, &
                          3,-5, 2], [3,3])) 
  b = [25,-10,-4]

  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
  x = solve(A,b) 

  print *, 'solution: ',x
  ! 5.0, 3.0, -2.0

end program example_solve1


program example_solve2
  use stdlib_linalg_constants, only: sp
  use stdlib_linalg, only: solve, linalg_state_type 
  implicit none

  complex(sp), allocatable :: A(:,:),b(:),x(:)

  ! Solve a system of 3 complex linear equations: 
  !  2x +      iy + 2z = (5-i)
  ! -ix + (4-3i)y + 6z = i
  !  4x +      3y +  z = 1

  ! Note: Fortran is column-major! -> transpose 
  A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), &
                         (0.0,-1.0),(4.0,-3.0),(6.0,0.0), &
                         (4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3])) 
  b = [(5.0,-1.0),(0.0,1.0),(1.0,0.0)]

  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
  x = solve(A,b) 

  print *, 'solution: ',x
  !   (1.0947,0.3674) (-1.519,-0.4539) (1.1784,-0.1078)

end program example_solve2

solve_lu - Solves a linear matrix equation or a linear system of equations (subroutine interface).

Status

Stable

Description

This subroutine computes the solution to a linear matrix equation $A \cdot x = b$, where $A$ is a square, full-rank, real or complex matrix.

Result vector or array x returns the exact solution to within numerical precision, provided that the matrix is not ill-conditioned. An error is returned if the matrix is rank-deficient or singular to working precision. If all optional arrays are provided by the user, no internal allocations take place. The solver is based on LAPACK's *GESV backends.

Syntax

Simple (Pure) interface:

call solve_lu (a, b, x)

Expert (Pure) interface:

call solve_lu (a, b, x [, pivot, overwrite_a, err])

Arguments

a: Shall be a rank-2 real or complex square array containing the coefficient matrix. It is normally an intent(in) argument. If overwrite_a=.true., it is an intent(inout) argument and is destroyed by the call.

b: Shall be a rank-1 or rank-2 array of the same kind as a, containing the right-hand-side vector(s). It is an intent(in) argument.

x: Shall be a rank-1 or rank-2 array of the same kind and size as b, that returns the solution(s) to the system. It is an intent(inout) argument, and must have the contiguous property.

pivot (optional): Shall be a rank-1 array of the same kind and matrix dimension as a, providing storage for the diagonal pivot indices. It is an intent(inout) arguments, and returns the diagonal pivot indices.

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.

Return value

For a full-rank matrix, returns an array value that represents the solution to the linear system of equations.

Raises LINALG_ERROR if the matrix is singular to working precision. Raises LINALG_VALUE_ERROR if the matrix and rhs vectors have invalid/incompatible sizes. If err is not present, exceptions trigger an error stop.

Example

program example_solve3
  use stdlib_linalg_constants, only: sp,ilp
  use stdlib_linalg, only: solve_lu, linalg_state_type 
  implicit none

  integer(ilp) :: test
  integer(ilp), allocatable :: pivot(:)
  complex(sp), allocatable :: A(:,:),b(:),x(:)

  ! Solve a system of 3 complex linear equations: 
  !  2x +      iy + 2z = (5-i)
  ! -ix + (4-3i)y + 6z = i
  !  4x +      3y +  z = 1

  ! Note: Fortran is column-major! -> transpose 
  A = transpose(reshape([(2.0, 0.0),(0.0, 1.0),(2.0,0.0), &
                         (0.0,-1.0),(4.0,-3.0),(6.0,0.0), &
                         (4.0, 0.0),(3.0, 0.0),(1.0,0.0)] , [3,3])) 

  ! Pre-allocate x
  allocate(b(size(A,2)),pivot(size(A,2)))
  allocate(x,mold=b)

  ! Call system many times avoiding reallocation 
  do test=1,100
     b = test*[(5.0,-1.0),(0.0,1.0),(1.0,0.0)]
     call solve_lu(A,b,x,pivot)
     print "(i3,'-th solution: ',*(1x,f12.6))", test,x
  end do

end program example_solve3

lstsq - Computes the least squares solution to a linear matrix equation.

Status

Stable

Description

This function computes the least-squares solution to a linear matrix equation $A \cdot x = b$.

Result vector x returns the approximate solution that minimizes the 2-norm $|| A \cdot x - b ||_2$, i.e., it contains the least-squares solution to the problem. Matrix A may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's *GELSD backends.

Syntax

x = lstsq (a, b, [, cond, overwrite_a, rank, err])

Arguments

a: Shall be a rank-2 real or complex array containing the coefficient matrix. It is an intent(inout) argument.

b: Shall be a rank-1 or rank-2 array of the same kind as a, containing one or more right-hand-side vector(s), each in its leading dimension. It is an intent(in) argument.

cond (optional): Shall be a scalar real value cut-off threshold for rank evaluation: s_i >= cond*maxval(s), i=1:rank. Shall be a scalar, intent(in) argument.

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.

rank (optional): Shall be an integer scalar value, that contains the rank of input matrix A. This is an intent(out) argument.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return value

Returns an array value of the same kind and rank as b, containing the solution(s) to the least squares system.

Raises LINALG_ERROR if the underlying Singular Value Decomposition process did not converge. Raises LINALG_VALUE_ERROR if the matrix and right-hand-side vector have invalid/incompatible sizes. Exceptions trigger an error stop.

Example

! Least-squares solver: functional interface
program example_lstsq1
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: lstsq
  implicit none

  integer, allocatable :: x(:),y(:)
  real(dp), allocatable :: A(:,:),b(:),coef(:)

  ! Data set
  x = [1, 2, 2]
  y = [5, 13, 25]

  ! Fit three points using a parabola, least squares method
  ! A = [1 x x**2]
  A = reshape([[1,1,1],x,x**2],[3,3])
  b = y

  ! Get coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
  coef = lstsq(A,b)

  print *, 'parabola: ',coef
  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 


end program example_lstsq1

solve_lstsq - Compute the least squares solution to a linear matrix equation (subroutine interface).

Status

Stable

Description

This subroutine computes the least-squares solution to a linear matrix equation $A \cdot x = b$.

Result vector x returns the approximate solution that minimizes the 2-norm $|| A \cdot x - b ||_2$, i.e., it contains the least-squares solution to the problem. Matrix A may be full-rank, over-determined, or under-determined. The solver is based on LAPACK's *GELSD backends.

Syntax

call solve_lstsq (a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])

Arguments

a: Shall be a rank-2 real or complex array containing the coefficient matrix. It is an intent(inout) argument.

b: Shall be a rank-1 or rank-2 array of the same kind as a, containing one or more right-hand-side vector(s), each in its leading dimension. It is an intent(in) argument.

x: Shall be an array of same kind and rank as b, and leading dimension of at least n, containing the solution(s) to the least squares system. It is an intent(inout) argument.

real_storage (optional): Shall be a real rank-1 array of the same kind a, providing working storage for the solver. It minimum size can be determined with a call to lstsq_space. It is an intent(inout) argument.

int_storage (optional): Shall be an integer rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to lstsq_space. It is an intent(inout) argument.

cmpl_storage (optional): For complex systems, it shall be a complex rank-1 array, providing working storage for the solver. It minimum size can be determined with a call to lstsq_space. It is an intent(inout) argument.

cond (optional): Shall be a scalar real value cut-off threshold for rank evaluation: s_i >= cond*maxval(s), i=1:rank. Shall be a scalar, intent(in) argument.

singvals (optional): Shall be a real rank-1 array of the same kind a and size at least min(m,n), returning the list of singular values s(i)>=cond*maxval(s) from the internal SVD, in descending order of magnitude. It is an intent(out) argument.

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.

rank (optional): Shall be an integer scalar value, that contains the rank of input matrix A. This is an intent(out) argument.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return value

Returns an array value that represents the solution to the least squares system.

Raises LINALG_ERROR if the underlying Singular Value Decomposition process did not converge. Raises LINALG_VALUE_ERROR if the matrix and right-hand-side vector have invalid/incompatible sizes. Exceptions trigger an error stop.

Example

! Demonstrate expert subroutine interface with pre-allocated arrays
program example_lstsq2
  use stdlib_linalg_constants, only: dp,ilp
  use stdlib_linalg, only: solve_lstsq, lstsq_space, linalg_state_type
  implicit none

  integer, allocatable :: x(:),y(:)
  real(dp), allocatable :: A(:,:),b(:),coef(:),real_space(:),singvals(:)
  integer(ilp), allocatable :: int_space(:)
  integer(ilp) :: lrwork,liwork,arank

  ! Data set
  x = [1, 2, 2]
  y = [5, 13, 25]

  ! Fit three points using a parabola, least squares method
  ! A = [1 x x**2]
  A = reshape([[1,1,1],x,x**2],[3,3])
  b = y

  ! Get storage sizes for the arrays and pre-allocate data
  call lstsq_space(A,b,lrwork,liwork)
  allocate(coef(size(x)),real_space(lrwork),int_space(liwork),singvals(minval(shape(A))))

  ! Solve coefficients of y = coef(1) + x*coef(2) + x^2*coef(3)
  ! with no internal allocations
  call solve_lstsq(A,b,x=coef,              &
                   real_storage=real_space, &
                   int_storage=int_space,   &
                   singvals=singvals,       &
                   overwrite_a=.true.,      &
                   rank=arank)

  print *, 'parabola: ',coef
  ! parabola:  -0.42857142857141695        1.1428571428571503        4.2857142857142811 
  print *, 'rank: ',arank
  ! rank: 2


end program example_lstsq2

lstsq_space - Compute internal working space requirements for the least squares solver.

Status

Stable

Description

This subroutine computes the internal working space requirements for the least-squares solver, solve_lstsq .

Syntax

call lstsq_space (a, b, lrwork, liwork [, lcwork])

Arguments

a: Shall be a rank-2 real or complex array containing the linear system coefficient matrix. It is an intent(in) argument.

b: Shall be a rank-1 or rank-2 array of the same kind as a, containing the system's right-hand-side vector(s). It is an intent(in) argument.

lrwork: Shall be an integer scalar, that returns the minimum array size required for the real working storage to this system.

liwork: Shall be an integer scalar, that returns the minimum array size required for the integer working storage to this system.

lcwork (complex a, b): For a complex system, shall be an integer scalar, that returns the minimum array size required for the complex working storage to this system.

det - Computes the determinant of a square matrix

Status

Stable

Description

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.

Syntax

c = det (a [, overwrite_a, err])

Arguments

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.

Return value

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.

Example

program example_determinant
  use stdlib_kinds, only: dp
  use stdlib_linalg, only: det, linalg_state_type
  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_determinant

.det. - Determinant operator of a square matrix

Status

Stable

Description

This operator returns the determinant of a real square matrix.

This interface is equivalent to the pure version of determinant det.

Syntax

c = [[stdlib_linalg(module):operator(.det.)(interface)]] a

Arguments

a: Shall be a rank-2 square array of any real or complex kinds. It is an intent(in) argument.

Return value

Returns a real scalar value 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 trigger an error stop.

Example

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

qr - Compute the QR factorization of a matrix

Status

Experimental

Description

This subroutine computes the QR factorization of a real or complex matrix: $A = Q R$ where $Q$ is orthonormal and $R$ is upper-triangular. Matrix $A$ has size [m,n], with $m \ge n$.

The results are returned in output matrices $Q$ and $R$, that have the same type and kind as $A$. Given k = min(m,n), one can write $A = \( Q_1 Q_2$ \cdot $\frac{R_1}{0}$ ). Because the lower rows of $R$ are zeros, a reduced problem $A = Q_1 R_1$ may be solved. The size of the input arguments determines what problem is solved: on full matrices (shape(Q)==[m,m], shape(R)==[m,n]), the full problem is solved. On reduced matrices (shape(Q)==[m,k], shape(R)==[k,n]), the reduced problem is solved.

Syntax

call qr (a, q, r, [, storage] [, overwrite_a] [, err])

Arguments

a: Shall be a rank-2 real or complex array containing the coefficient matrix of size [m,n]. It is an intent(in) argument, if overwrite_a=.false.. Otherwise, it is an intent(inout) argument, and is destroyed upon return.

q: Shall be a rank-2 array of the same kind as a, containing the orthonormal matrix q. It is an intent(out) argument. It should have a shape equal to either [m,m] or [m,k], whether the full or the reduced problem is sought for.

r: Shall be a rank-2 array of the same kind as a, containing the upper triangular matrix r. It is an intent(out) argument. It should have a shape equal to either [m,n] or [k,n], whether the full or the reduced problem is sought for.

storage (optional): Shall be a rank-1 array of the same type and kind as a, providing working storage for the solver. Its minimum size can be determined with a call to qr_space. It is an intent(out) argument.

overwrite_a (optional): Shall be an input logical flag (default: .false.). If .true., input matrix a will be used as temporary storage and overwritten. This avoids internal data allocation. It is an intent(in) argument.

err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument.

Return value

Returns the QR factorization matrices into the $Q$ and $R$ arguments.

Raises LINALG_VALUE_ERROR if any of the matrices has invalid or unsuitable size for the full/reduced problem. Raises LINALG_ERROR on insufficient user storage space. If the state argument err is not present, exceptions trigger an error stop.

Example

program example_qr  
  use stdlib_linalg, only: qr
  implicit none(type,external)
  real :: A(104, 32), Q(104,32), R(32,32)

  ! Create a random matrix
  call random_number(A)

  ! Compute its QR factorization (reduced)
  call qr(A,Q,R)

  ! Test factorization: Q*R = A 
  print *, maxval(abs(matmul(Q,R)-A)) 

end program example_qr

qr_space - Compute internal working space requirements for the QR factorization.

Status

Experimental

Description

This subroutine computes the internal working space requirements for the QR factorization, qr .

Syntax

call qr_space (a, lwork, [, err])

Arguments

a: Shall be a rank-2 real or complex array containing the coefficient matrix. It is an intent(in) argument.

lwork: Shall be an integer scalar, that returns the minimum array size required for the working storage in qr to factorize a.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Example

! QR example with pre-allocated storage
program example_qr_space
  use stdlib_linalg_constants, only: ilp
  use stdlib_linalg, only: qr, qr_space, linalg_state_type
  implicit none(type,external)
  real :: A(104, 32), Q(104,32), R(32,32)
  real, allocatable :: work(:)
  integer(ilp) :: lwork
  type(linalg_state_type) :: err

  ! Create a random matrix
  call random_number(A)

  ! Prepare QR workspace
  call qr_space(A,lwork)
  allocate(work(lwork))

  ! Compute its QR factorization (reduced)
  call qr(A,Q,R,storage=work,err=err)

  ! Test factorization: Q*R = A 
  print *, maxval(abs(matmul(Q,R)-A))
  print *, err%print() 

end program example_qr_space

eig - Eigenvalues and Eigenvectors of a Square Matrix

Status

Stable

Description

This subroutine computes the solution to the eigenproblem $A \cdot \bar{v} - \lambda \cdot \bar{v}$, where $A$ is a square, full-rank, real or complex matrix.

Result array lambda returns the eigenvalues of $A$. The user can request eigenvectors to be returned: if provided, on output left will contain the left eigenvectors, right the right eigenvectors of $A$. Both left and right are rank-2 arrays, where eigenvectors are stored as columns. The solver is based on LAPACK's *GEEV backends.

Syntax

call eig (a, lambda [, right] [,left] [,overwrite_a] [,err])

Arguments

a : real or complex square array containing the coefficient matrix. If overwrite_a=.false., it is an intent(in) argument. Otherwise, it is an intent(inout) argument and is destroyed by the call.

lambda: Shall be a complex or real rank-1 array of the same kind as a, containing the eigenvalues, or their real component only. It is an intent(out) argument.

right (optional): Shall be a complex rank-2 array of the same size and kind as a, containing the right eigenvectors of a. It is an intent(out) argument.

left (optional): Shall be a complex rank-2 array of the same size and kind as a, containing the left eigenvectors of a. It is an intent(out) argument.

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.

Return value

Raises LINALG_ERROR if the calculation did not converge. Raises LINALG_VALUE_ERROR if any matrix or arrays have invalid/incompatible sizes. Raises LINALG_VALUE_ERROR if the real component is only requested, but the eigenvalues have non-trivial imaginary parts. If err is not present, exceptions trigger an error stop.

Example

! Eigendecomposition of a real square matrix 
program example_eig
  use stdlib_linalg, only: eig
  implicit none

  integer :: i
  real, allocatable :: A(:,:)
  complex, allocatable :: lambda(:),vectors(:,:)

  ! Decomposition of this square matrix
  ! NB Fortran is column-major -> transpose input
  A = transpose(reshape( [ [2, 2, 4], &
                           [1, 3, 5], &
                           [2, 3, 4] ], [3,3] )) 

  ! Get eigenvalues and right eigenvectors
  allocate(lambda(3),vectors(3,3))

  call eig(A, lambda, right=vectors)

  do i=1,3
     print *, 'eigenvalue  ',i,': ',lambda(i)
     print *, 'eigenvector ',i,': ',vectors(:,i)
  end do

end program example_eig

eigh - Eigenvalues and Eigenvectors of a Real symmetric or Complex Hermitian Square Matrix

Status

Stable

Description

This subroutine computes the solution to the eigendecomposition $A \cdot \bar{v} - \lambda \cdot \bar{v}$, where $A$ is a square, full-rank, real symmetric $A = A^T$ or complex Hermitian $A = A^H$ matrix.

Result array lambda returns the real eigenvalues of $A$. The user can request the orthogonal eigenvectors to be returned: on output vectors may contain the matrix of eigenvectors, returned as a column.

Normally, only the lower triangular part of $A$ is accessed. On input, logical flag upper_a allows the user to request what triangular part of the matrix should be used.

The solver is based on LAPACK's *SYEV and *HEEV backends.

Syntax

call eigh (a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])

Arguments

a : real or complex square array containing the coefficient matrix. It is an intent(in) argument. If overwrite_a=.true., it is an intent(inout) argument and is destroyed by the call.

lambda: Shall be a complex rank-1 array of the same precision as a, containing the eigenvalues. It is an intent(out) argument.

vectors (optional): Shall be a rank-2 array of the same type, size and kind as a, containing the eigenvectors of a. It is an intent(out) argument.

upper_a (optional): Shall be an input logical flag. If .true., the upper triangular part of a will be accessed. Otherwise, the lower triangular part will be accessed. It is an intent(in) argument.

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.

Return value

Raises LINALG_ERROR if the calculation did not converge. Raises LINALG_VALUE_ERROR if any matrix or arrays have invalid/incompatible sizes. If err is not present, exceptions trigger an error stop.

Example

! Eigendecomposition of a real symmetric matrix 
program example_eigh
  use stdlib_linalg, only: eigh
  implicit none

  integer :: i
  real, allocatable :: A(:,:),lambda(:),vectors(:,:)
  complex, allocatable :: cA(:,:),cvectors(:,:)

  ! Decomposition of this symmetric matrix
  ! NB Fortran is column-major -> transpose input
  A = transpose(reshape( [ [2, 1, 4], &
                           [1, 3, 5], &
                           [4, 5, 4] ], [3,3] )) 

  ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors
  allocate(lambda(3),vectors(3,3))  
  call eigh(A, lambda, vectors=vectors)

  print *, 'Real matrix'
  do i=1,3
     print *, 'eigenvalue  ',i,': ',lambda(i)
     print *, 'eigenvector ',i,': ',vectors(:,i)
  end do

  ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors
  cA = A

  allocate(cvectors(3,3))  
  call eigh(cA, lambda, vectors=cvectors)

  print *, 'Complex matrix'
  do i=1,3
     print *, 'eigenvalue  ',i,': ',lambda(i)
     print *, 'eigenvector ',i,': ',cvectors(:,i)
  end do  

end program example_eigh

eigvals - Eigenvalues of a Square Matrix

Status

Stable

Description

This function returns the eigenvalues to matrix $A$: a square, full-rank, real or complex matrix. The eigenvalues are solutions to the eigenproblem $A \cdot \bar{v} - \lambda \cdot \bar{v}$.

Result array lambda is complex, and returns the eigenvalues of $A$. The solver is based on LAPACK's *GEEV backends.

Syntax

lambda = eigvals (a, [,err])

Arguments

a : real or complex square array containing the coefficient matrix. It is an intent(in) argument.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return value

Returns a complex array containing the eigenvalues of a.

Raises LINALG_ERROR if the calculation did not converge. Raises LINALG_VALUE_ERROR if any matrix or arrays have invalid/incompatible sizes. If err is not present, exceptions trigger an error stop.

Example

! Eigenvalues of a general real / complex matrix 
program example_eigvals
  use stdlib_linalg, only: eigvals
  implicit none

  real, allocatable :: A(:,:),lambda(:)
  complex, allocatable :: cA(:,:),clambda(:)

  ! NB Fortran is column-major -> transpose input
  A = transpose(reshape( [ [2, 8, 4], &
                           [1, 3, 5], &
                           [9, 5,-2] ], [3,3] )) 

  ! Note: real symmetric matrix
  lambda = eigvals(A)
  print *, 'Real    matrix eigenvalues: ',lambda

  ! Complex general matrix
  cA = cmplx(A, -2*A)
  clambda = eigvals(cA)
  print *, 'Complex matrix eigenvalues: ',clambda

end program example_eigvals

eigvalsh - Eigenvalues of a Real Symmetric or Complex Hermitian Square Matrix

Status

Stable

Description

This function returns the eigenvalues to matrix $A$: a where $A$ is a square, full-rank, real symmetric $A = A^T$ or complex Hermitian $A = A^H$ matrix. The eigenvalues are solutions to the eigenproblem $A \cdot \bar{v} - \lambda \cdot \bar{v}$.

Result array lambda is real, and returns the eigenvalues of $A$. The solver is based on LAPACK's *SYEV and *HEEV backends.

Syntax

lambda = eigvalsh (a, [, upper_a] [,err])

Arguments

a : real or complex square array containing the coefficient matrix. It is an intent(in) argument.

upper_a (optional): Shall be an input logical flag. If .true., the upper triangular part of a will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an intent(in) argument.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return value

Returns a real array containing the eigenvalues of a.

Raises LINALG_ERROR if the calculation did not converge. Raises LINALG_VALUE_ERROR if any matrix or arrays have invalid/incompatible sizes. If err is not present, exceptions trigger an error stop.

Example

! Eigenvalues of a real symmetric / complex hermitian matrix 
program example_eigvalsh
  use stdlib_linalg, only: eigvalsh
  implicit none

  real, allocatable :: A(:,:),lambda(:)
  complex, allocatable :: cA(:,:)

  ! Decomposition of this symmetric matrix
  ! NB Fortran is column-major -> transpose input
  A = transpose(reshape( [ [2, 1, 4], &
                           [1, 3, 5], &
                           [4, 5, 4] ], [3,3] )) 

  ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors
  lambda = eigvalsh(A)
  print *, 'Symmetric matrix eigenvalues: ',lambda

  ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors
  cA = A
  lambda = eigvalsh(cA)
  print *, 'Hermitian matrix eigenvalues: ',lambda

end program example_eigvalsh

svd - Compute the singular value decomposition of a rank-2 array (matrix).

Status

Stable

Description

This subroutine computes the singular value decomposition of a real or complex rank-2 array (matrix) $A = U \cdot S \cdot \V^T$. The solver is based on LAPACK's *GESDD backends.

Result vector s returns the array of singular values on the diagonal of $S$. If requested, u contains the left singular vectors, as columns of $U$. If requested, vt contains the right singular vectors, as rows of $V^T$.

Syntax

call svd (a, s, [, u, vt, overwrite_a, full_matrices, err])

Class

Subroutine

Arguments

a: Shall be a rank-2 real or complex array containing the coefficient matrix of size [m,n]. It is an intent(inout) argument, but returns unchanged unless overwrite_a=.true..

s: Shall be a rank-1 real array, returning the list of k = min(m,n) singular values. It is an intent(out) argument.

u (optional): Shall be a rank-2 array of same kind as a, returning the left singular vectors of a as columns. Its size should be [m,m] unless full_matrices=.false., in which case, it can be [m,min(m,n)]. It is an intent(out) argument.

vt (optional): Shall be a rank-2 array of same kind as a, returning the right singular vectors of a as rows. Its size should be [n,n] unless full_matrices=.false., in which case, it can be [min(m,n),n]. It is an intent(out) argument.

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. By default, overwrite_a=.false.. It is an intent(in) argument.

full_matrices (optional): Shall be an input logical flag. If .true. (default), matrices u and vt shall be full-sized. Otherwise, their secondary dimension can be resized to min(m,n). See u, v for details.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return values

Returns an array s that contains the list of singular values of matrix a. If requested, returns a rank-2 array u that contains the left singular vectors of a along its columns. If requested, returns a rank-2 array vt that contains the right singular vectors of a along its rows.

Raises LINALG_ERROR if the underlying Singular Value Decomposition process did not converge. Raises LINALG_VALUE_ERROR if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an error stop, unless argument err is present.

Example

svdvals - Compute the singular values of a rank-2 array (matrix).

Status

Stable

Description

This subroutine computes the singular values of a real or complex rank-2 array (matrix) from its singular value decomposition $A = U \cdot S \cdot \V^T$. The solver is based on LAPACK's *GESDD backends.

Result vector s returns the array of singular values on the diagonal of $S$.

Syntax

s = svdvals (a [, err])

Arguments

a: Shall be a rank-2 real or complex array containing the coefficient matrix of size [m,n]. It is an intent(in) argument.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return values

Returns an array s that contains the list of singular values of matrix a.

Raises LINALG_ERROR if the underlying Singular Value Decomposition process did not converge. Raises LINALG_VALUE_ERROR if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an error stop, unless argument err is present.

Example

! Singular Values
program example_svdvals
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: svdvals
  implicit none

  real(dp), allocatable :: A(:,:),s(:)
  character(*), parameter :: fmt="(a,*(1x,f12.8))"

  ! We want to find the singular values of matrix: 
  ! 
  !    A = [ 3  2  2]
  !        [ 2  3 -2]  
  !   
  A = transpose(reshape([ 3, 2, 2, &
                          2, 3,-2], [3,2]))

  ! Get singular values
  s = svdvals(A)

  ! Singular values: [5, 3]
  print fmt, '    '
  print fmt, 'S = ',s
  print fmt, '    '

end program example_svdvals

cholesky - Compute the Cholesky factorization of a rank-2 square array (matrix)

Status

Experimental

Description

This subroutine computes the Cholesky factorization of a real or complex rank-2 square array (matrix), $A = L \cdot L^T$, or $A = U^T \cdot U$. $A$ is symmetric or complex Hermitian, and $L$, $U$ are lower- or upper-triangular, respectively. The solver is based on LAPACK's *POTRF backends.

Syntax

call cholesky (a, c, lower [, other_zeroed] [, err])

Class

Subroutine

Arguments

a: Shall be a rank-2 square real or complex array containing the coefficient matrix of size [n,n]. It is an intent(inout) argument, but returns unchanged if the argument c is present.

c (optional): Shall be a rank-2 square real or complex of the same size and kind as a. It is an intent(out) argument, that returns the triangular Cholesky matrix L or U.

lower: Shall be an input logical flag. If .true., the lower triangular decomposition $A = L \cdot L^T$ will be performed. If .false., the upper decomposition $A = U^T \cdot U$ will be performed.

other_zeroed (optional): Shall be an input logical flag. If .true., the unused part of the output matrix will contain zeroes. Otherwise, it will not be accessed. This saves cpu time. By default, other_zeroed=.true.. It is an intent(in) argument.

err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument.

Return values

The factorized matrix is returned in-place overwriting a if no other arguments are provided. Otherwise, it can be provided as a second argument c. In this case, a is not overwritten. The logical flag lower determines whether the lower- or the upper-triangular factorization should be performed.

Results are returned on the applicable triangular region of the output matrix, while the unused triangular region is filled by zeroes by default. Optional argument other_zeroed, if .false. allows the expert user to avoid zeroing the unused part; however, in this case, the unused region of the matrix is not accessed and will usually contain invalid values.

Raises LINALG_ERROR if the underlying process did not converge. Raises LINALG_VALUE_ERROR if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an error stop, unless argument err is present.

Example

! Cholesky factorization: subroutine interface 
program example_cholesky
  use stdlib_linalg, only: cholesky
  implicit none

  real, dimension(3,3) :: A,L,U 

  ! Set real matrix
  A = reshape( [ [6, 15, 55], &
                 [15, 55, 225], &
                 [55, 225, 979] ], [3,3] )

  ! Decompose (lower)
  call cholesky(A, L, lower=.true.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(L,transpose(L))))

  ! Decompose (upper)
  call cholesky(A, U, lower=.false.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(transpose(U),U)))

end program example_cholesky

chol - Compute the Cholesky factorization of a rank-2 square array (matrix)

Status

Experimental

Description

This is a pure functional interface for the Cholesky factorization of a real or complex rank-2 square array (matrix) computed as $A = L \cdot L^T$, or $A = U^T \cdot U$. $A$ is symmetric or complex Hermitian, and $L$, $U$ are lower- or upper-triangular, respectively. The solver is based on LAPACK's *POTRF backends.

Result matrix c has the same size and kind as a, and returns the lower or upper triangular factor.

Syntax

c = chol (a, lower [, other_zeroed])

Arguments

a: Shall be a rank-2 square real or complex array containing the coefficient matrix of size [n,n]. It is an intent(inout) argument, but returns unchanged if argument c is present.

lower: Shall be an input logical flag. If .true., the lower triangular decomposition $A = L \cdot L^T$ will be performed. If .false., the upper decomposition $A = U^T \cdot U$ will be performed.

other_zeroed (optional): Shall be an input logical flag. If .true., the unused part of the output matrix will contain zeroes. Otherwise, it will not be accessed. This saves cpu time. By default, other_zeroed=.true.. It is an intent(in) argument.

Return values

Returns a rank-2 array c of the same size and kind as a, that contains the triangular Cholesky matrix L or U.

Raises LINALG_ERROR if the underlying process did not converge. Raises LINALG_VALUE_ERROR if the matrix or any of the output arrays invalid/incompatible sizes. Exceptions trigger an error stop, unless argument err is present.

Example

! Cholesky factorization: function interface 
program example_chol
  use stdlib_linalg, only: chol
  implicit none

  real, allocatable, dimension(:,:) :: A,L,U 

  ! Set real matrix
  A = reshape( [ [6, 15, 55], &
                 [15, 55, 225], &
                 [55, 225, 979] ], [3,3] )

  ! Decompose (lower)
  L = chol(A, lower=.true.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(L,transpose(L))))

  ! Decompose (upper)
  U = chol(A, lower=.false.)

  ! Compare decomposition 
  print *, maxval(abs(A-matmul(transpose(U),U)))

end program example_chol

.inv. - Inverse operator of a square matrix

Status

Stable

Description

This operator returns the inverse of a real or complex square matrix $A$. The inverse $A^{-1}$ is defined such that $A \cdot A^{-1} = A^{-1} \cdot A = I_n$.

This interface is equivalent to the function inv.

Syntax

b = [[stdlib_linalg(module):operator(.inv.)(interface)]] a

Arguments

a: Shall be a rank-2 square array of any real or complex kinds. It is an intent(in) argument.

Return value

Returns a rank-2 square array with the same type, kind and rank as a, that contains the inverse of a.

If an exception occurred on input errors, or singular matrix, NaNs will be returned. For fine-grained error control in case of singular matrices prefer the subroutine and the function interfaces.

Example

! Matrix inversion example: operator interface
program example_inverse_operator
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: operator(.inv.),eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))

  ! Invert matrix
  Am1 = .inv.A

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_operator

invert - Inversion of a square matrix

Status

Stable

Description

This subroutine inverts a square real or complex matrix in-place. The inverse $A^{-1}$ is defined such that $A \cdot A^{-1} = A^{-1} \cdot A = I_n$.

On return, the input matrix a is replaced by its inverse. The solver is based on LAPACK's *GETRF and *GETRI backends.

Syntax

call invert (a, [,inva] [, pivot] [, err])

Arguments

a: Shall be a rank-2, square, real or complex array containing the coefficient matrix. If inva is provided, it is an intent(in) argument. If inva is not provided, it is an intent(inout) argument: on output, it is replaced by the inverse of a.

inva (optional): Shall be a rank-2, square, real or complex array with the same size, and kind as a. On output, it contains the inverse of a.

pivot (optional): Shall be a rank-1 array of the same kind and matrix dimension as a, that contains the diagonal pivot indices on return. It is an intent(inout) argument.

err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.

Return value

Computes the inverse of the matrix $A$, $A^{-1}\, and returns it either in \( A$ or in another matrix.

Raises LINALG_ERROR if the matrix is singular or has invalid size. Raises LINALG_VALUE_ERROR if inva and a do not have the same size. If err is not present, exceptions trigger an error stop.

Example

! Matrix inversion example: in-place inversion
program example_inverse_inplace
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: invert,eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))
  Am1 = A

  ! Invert matrix
  call invert(Am1)

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_inplace

! Matrix inversion example: subroutine interface
program example_inverse_subroutine
  use stdlib_linalg_constants, only: dp
  use stdlib_linalg, only: invert,eye
  implicit none

  real(dp) :: A(2,2), Am1(2,2)

  ! Input matrix (NB Fortran is column major! input columns then transpose)
  A = transpose(reshape( [4, 3, &
                          3, 2], [2,2] ))

  ! Invert matrix
  call invert(A,Am1)

  print *, '         |',Am1(1,:),'|' ! | -2  3 |
  print *, ' inv(A)= |',Am1(2,:),'|' ! |  3 -4 |

  ! Final check 
  print *, 'CHECK passed? ',matmul(A,Am1)==eye(2)

end program example_inverse_subroutine