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
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).
program example_gemv
use stdlib_linalg, only: eye
use stdlib_linalg_blas, only: sp,gemv
implicit none
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
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 arrayStable
Pure function.
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 matrixStable
Pure function.
Constructs the identity matrix.
I = stdlib_linalg (dim1 [, dim2] [, mold])
dim1: A scalar of type integer. This is an intent(in) argument and specifies the number of rows.dim2: A scalar of type integer. This is an optional intent(in) argument specifying the number of columns. If not provided, the matrix is square (dim1 = dim2).mold: A scalar of any supported integer, real, or complex type. This is an optional intent(in) argument. If provided, the returned identity matrix will have the same type and kind as mold. If not provided, the matrix will be of type real(real64) by default.Returns the identity matrix, with ones on the main diagonal and zeros elsewhere.
real(real64), which is recommended for arithmetic safety.mold argument is provided, the return value will match the type and kind of mold, allowing for arbitrary integer, real, or complex return types.!> Return default type (real64)
A = eye(2,2)/2 !! A == diag([0.5_dp, 0.5_dp])
!> Return 32-bit complex
A = eye(2,2, mold=(0.0,0.0))/2 !! A == diag([(0.5,0.5), (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 matrixStable
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
hermitian - Compute the Hermitian version of a rank-2 matrixExperimental
Compute the Hermitian version of a rank-2 matrix.
For complex matrices, the function returns the conjugate transpose (conjg(transpose(a))).
For real or integer matrices, the function returns the transpose (transpose(a)).
h = hermitian (a)
a: Shall be a rank-2 array of type integer, real, or complex. The input matrix a is not modified.
Returns a rank-2 array of the same shape and type as a. If a is of type complex, the Hermitian matrix is computed as conjg(transpose(a)).
For real or integer types, it is equivalent to the intrinsic transpose(a).
! Example program demonstrating the usage of hermitian
program example_hermitian
use stdlib_linalg, only: hermitian
implicit none
complex, dimension(2, 2) :: A, AT
! Define input matrices
A = reshape([(1,2),(3,4),(5,6),(7,8)],[2,2])
! Compute Hermitian matrices
AT = hermitian(A)
print *, "Original Complex Matrix:"
print *, A
print *, "Hermitian Complex Matrix:"
print *, AT
end program example_hermitian
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
solve - Solves a linear matrix equation or a linear system of equations.Stable
This function computes the solution to a linear matrix equation , where 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.
Pure interface:
x = solve (a, b)
Expert interface:
x = solve (a, b [, overwrite_a], err)
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.
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.
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).Stable
This subroutine computes the solution to a linear matrix equation , where 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.
Simple (Pure) interface:
call solve_lu (a, b, x)
Expert (Pure) interface:
call solve_lu (a, b, x [, pivot, overwrite_a, err])
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.
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.
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.Stable
This function computes the least-squares solution to a linear matrix equation .
Result vector x returns the approximate solution that minimizes the 2-norm , 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.
x = lstsq (a, b, [, cond, overwrite_a, rank, err])
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.
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.
! 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).Stable
This subroutine computes the least-squares solution to a linear matrix equation .
Result vector x returns the approximate solution that minimizes the 2-norm , 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.
call solve_lstsq (a, b, x, [, real_storage, int_storage, [cmpl_storage, ] cond, singvals, overwrite_a, rank, err])
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.
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.
! 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.Stable
This subroutine computes the internal working space requirements for the least-squares solver, solve_lstsq .
call lstsq_space (a, b, lrwork, liwork [, lcwork])
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 matrixStable
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
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 matrixStable
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 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.
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 matrixExperimental
This subroutine computes the QR factorization of a real or complex matrix: where
is orthonormal and is upper-triangular. Matrix has size [m,n], with .
The results are returned in output matrices and , that have the same type and kind as .
Given k = min(m,n), one can write \cdot ).
Because the lower rows of are zeros, a reduced problem 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.
call qr (a, q, r, [, storage] [, overwrite_a] [, err])
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.
Returns the QR factorization matrices into the and 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.
program example_qr
use stdlib_linalg, only: qr
implicit none
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.Experimental
This subroutine computes the internal working space requirements for the QR factorization, qr .
call qr_space (a, lwork, [, err])
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.
! 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
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
schur - Compute the Schur decomposition of a matrixExperimental
This subroutine computes the Schur decomposition of a real or complex matrix: , where is unitary (orthonormal) and is upper-triangular (for complex matrices) or quasi-upper-triangular (for real matrices, with possible blocks on the diagonal). Matrix has size [n,n].
The results are returned in output matrices and . Matrix is the Schur form, and matrix is the unitary transformation matrix such that . If requested, the eigenvalues of can also be returned as a complex array of size [n].
call schur (a, t [, z,] [, eigvals] [, overwrite_a] [, storage] [, err])
a: Shall be a rank-2 real or complex array containing the matrix to be decomposed. It is an intent(inout) argument if overwrite_a = .true.; otherwise, it is an intent(in) argument.
t: Shall be a rank-2 array of the same kind as a, containing the Schur form of the matrix. It is an intent(out) argument and should have a shape of [n,n].
z: Shall be a rank-2 array of the same kind as a, containing the unitary matrix . It is an intent(out) argument and is optional. If provided, it should have the shape [n,n].
eigvals (optional): Shall be a rank-1 complex or real array of the same kind as a, containing the eigenvalues of (the diagonal elements of ), or their real component only. The array must be of size [n]. If not provided, the eigenvalues are not returned. It is an intent(out) argument.
overwrite_a (optional): Shall be a logical flag (default: .false.). If .true., the input matrix a will be overwritten and destroyed upon return, avoiding internal data allocation. It is an intent(in) argument.
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 schur_space. It is an intent(inout) argument.
err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument. If not provided, exceptions trigger an error stop.
Returns the Schur decomposition matrices into the and arguments. If eigvals is provided, it will also return the eigenvalues of the matrix .
Raises LINALG_VALUE_ERROR if any of the matrices have invalid or unsuitable size for the decomposition.
Raises LINALG_VALUE_ERROR if the real component is only requested, but the eigenvalues have non-trivial imaginary parts.
Raises LINALG_ERROR on insufficient user storage space.
If the state argument err is not present, exceptions trigger an error stop.
! This example includes eigenvalue computation in addition to
! the Schur decomposition for a randomly generated matrix.
program example_schur_eigenvalues
use stdlib_linalg, only: schur
use stdlib_linalg_constants, only: dp
implicit none
integer, parameter :: n = 5
real(dp), dimension(n,n) :: A, T, Z
complex(dp), dimension(n) :: eigenvalues
! Create a random real-valued square matrix
call random_number(A)
! Compute the Schur decomposition and eigenvalues
call schur(A, T, Z, eigenvalues)
! Output results
print *, "Random Matrix A:"
print *, A
print *, "Schur Form Matrix T:"
print *, T
print *, "Orthogonal Matrix Z:"
print *, Z
print *, "Eigenvalues:"
print *, eigenvalues
! Test factorization: Z*T*Z^T = A
print *, "Max error in reconstruction:", maxval(abs(matmul(Z, matmul(T, transpose(Z))) - A))
end program example_schur_eigenvalues
schur_space - Compute internal working space requirements for the Schur decompositionExperimental
This subroutine computes the internal working space requirements for the Schur decomposition, schur.
call schur_space (a, lwork, [, err])
a: Shall be a rank-2 real or complex array containing the matrix to be decomposed. It is an intent(in) argument.
lwork: Shall be an integer scalar that returns the minimum array size required for the working storage in schur to decompose a. It is an intent(out) argument.
err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument.
Returns the required working storage size lwork for the Schur decomposition. This value can be used to pre-allocate a workspace array in case multiple Schur decompositions of the same matrix size are needed. If pre-allocated working arrays are provided, no internal memory allocations will take place during the decomposition.
eig - Eigenvalues and Eigenvectors of a Square MatrixStable
This subroutine computes the solution to the eigenproblem ,
where is a square, full-rank, real or complex matrix, or to the generalized eigenproblem ,
where is a square matrix with the same type, kind and size as .
Result array lambda returns the eigenvalues of . The user can request eigenvectors to be returned: if provided, on output left will contain the left eigenvectors, right the right eigenvectors of .
Both left and right are rank-2 arrays, where eigenvectors are stored as columns.
The solver is based on LAPACK's *GEEV (standard eigenproblem) and *GGEV (generalized eigenproblem) backends.
For the standard eigenproblem:
call eig (a, lambda [, right] [,left] [,overwrite_a] [,err])
For the generalized eigenproblem:
call eig `(a, b, lambda [, right] [, left] [, overwrite_a] [, overwrite_b] [, err])
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.
b: real or complex square array containing the second coefficient matrix. If overwrite_b=.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.
overwrite_b (optional): Shall be an input logical flag. If .true., input matrix b 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.
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.
! 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 MatrixStable
This subroutine computes the solution to the eigendecomposition ,
where is a square, full-rank, real symmetric or complex Hermitian matrix.
Result array lambda returns the real eigenvalues of . 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 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.
call eigh (a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])
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.
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.
! 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 MatrixStable
This function computes the eigenvalues for either a standard or generalized eigenproblem:
real or complex matrix.The eigenvalues are stored in the result array lambda, which is complex (even for real input matrices).
The solver uses LAPACK's *GEEV and *GGEV backends for the standard and generalized problems, respectively.
For the standard eigenproblem:
lambda = eigvals (a [, err])
For the generalized eigenproblem:
lambda = eigvals (a, b [, err])
a:
Shall be a real or complex square array containing the coefficient matrix. It is an intent(in) argument.
b (optional):
Shall be a real or complex square array containing the second coefficient matrix for the generalized problem. It is an intent(in) argument.
err (optional):
Shall be a type(linalg_state_type) value. This is an intent(out) argument.
Returns a complex rank-1 array containing the eigenvalues of the problem.
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.
! 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 MatrixStable
This function returns the eigenvalues to matrix : a where is a square, full-rank,
real symmetric or complex Hermitian matrix.
The eigenvalues are solutions to the eigenproblem .
Result array lambda is real, and returns the eigenvalues of .
The solver is based on LAPACK's *SYEV and *HEEV backends.
lambda = eigvalsh (a, [, upper_a] [,err])
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.
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.
! 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).Stable
This subroutine computes the singular value decomposition of a real or complex rank-2 array (matrix) .
The solver is based on LAPACK's *GESDD backends.
Result vector s returns the array of singular values on the diagonal of .
If requested, u contains the left singular vectors, as columns of .
If requested, vt contains the right singular vectors, as rows of .
call svd (a, s, [, u, vt, overwrite_a, full_matrices, err])
Subroutine
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.
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.
svdvals - Compute the singular values of a rank-2 array (matrix).Stable
This subroutine computes the singular values of a real or complex rank-2 array (matrix) from its singular
value decomposition . The solver is based on LAPACK's *GESDD backends.
Result vector s returns the array of singular values on the diagonal of .
s = svdvals (a [, err])
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.
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.
! 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)Experimental
This subroutine computes the Cholesky factorization of a real or complex rank-2 square array (matrix),
, or . is symmetric or complex Hermitian, and ,
are lower- or upper-triangular, respectively.
The solver is based on LAPACK's *POTRF backends.
call cholesky (a, c, lower [, other_zeroed] [, err])
Subroutine
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 will be performed. If .false., the upper decomposition 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.
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.
! 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)Experimental
This is a pure functional interface for the Cholesky factorization of a real or
complex rank-2 square array (matrix) computed as , or .
is symmetric or complex Hermitian, and , 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.
c = chol (a, lower [, other_zeroed])
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 will be performed. If .false., the upper decomposition 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.
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.
! 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 matrixStable
This operator returns the inverse of a real or complex square matrix .
The inverse is defined such that .
This interface is equivalent to the function inv.
b = [[stdlib_linalg(module):operator(.inv.)(interface)]] a
a: Shall be a rank-2 square array of any real or complex kinds. It is an intent(in) argument.
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.
! 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 matrixStable
This subroutine inverts a square real or complex matrix in-place.
The inverse is defined such that .
On return, the input matrix a is replaced by its inverse.
The solver is based on LAPACK's *GETRF and *GETRI backends.
call invert (a, [,inva] [, pivot] [, err])
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.
Computes the inverse of the matrix , 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.
! 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
inv - Inverse of a square matrix.Stable
This function returns the inverse of a square real or complex matrix in-place.
The inverse, , is defined such that .
The solver is based on LAPACK's *GETRF and *GETRI backends.
b inv (a, [, err])
a: Shall be a rank-2, square, real or complex array containing the coefficient matrix. It is an intent(inout) argument.
err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument.
Returns an array value of the same type, kind and rank as a, that contains the inverse matrix .
Raises LINALG_ERROR if the matrix is singular or has invalid size.
If err is not present, exceptions trigger an error stop.
! Matrix inversion example: function interface
program example_inverse_function
use stdlib_linalg_constants, only: dp
use stdlib_linalg, only: 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_function
pinv - Moore-Penrose pseudo-inverse of a matrixExperimental
This function computes the Moore-Penrose pseudo-inverse of a real or complex matrix.
The pseudo-inverse, , generalizes the matrix inverse and satisfies the conditions:
-
-
-
-
The computation is based on singular value decomposition (SVD). Singular values below a relative tolerance threshold , where is the largest singular value, are treated as zero.
b = pinv (a, [, rtol, err])
a: Shall be a rank-2, real or complex array of shape [m, n] containing the coefficient matrix.
It is an intent(in) argument.
rtol (optional): Shall be a scalar real value specifying the relative tolerance for singular value cutoff.
If rtol is not provided, the default relative tolerance is ,
where is the machine precision for the element type of a. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument.
Returns an array value of the same type, kind, and rank as a with shape [n, m], that contains the pseudo-inverse matrix .
Raises LINALG_ERROR if the underlying SVD did not converge.
Raises LINALG_VALUE_ERROR if a has invalid size.
If err is not present, exceptions trigger an error stop.
! Matrix pseudo-inversion example: function, subroutine, and operator interfaces
program example_pseudoinverse
use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type
implicit none
real :: A(15,5), Am1(5,15)
type(linalg_state_type) :: state
! Generate random matrix A (15x15)
call random_number(A)
! Pseudo-inverse: Function interfcae
Am1 = pinv(A, err=state)
print *, 'Max error (function) : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! User threshold
Am1 = pinv(A, rtol=0.001, err=state)
print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! Pseudo-inverse: Subroutine interface
call pseudoinvert(A, Am1, err=state)
print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! Operator interface
Am1 = .pinv.A
print *, 'Max error (operator) : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
end program example_pseudoinverse
pseudoinvert - Moore-Penrose pseudo-inverse of a matrixExperimental
This subroutine computes the Moore-Penrose pseudo-inverse of a real or complex matrix.
The pseudo-inverse is a generalization of the matrix inverse and satisfies the following properties:
-
-
-
-
The computation is based on singular value decomposition (SVD). Singular values below a relative tolerance threshold , where is the largest singular value, are treated as zero.
On return, matrix pinva [n, m] will store the pseudo-inverse of a [m, n].
call pseudoinvert (a, pinva [, rtol] [, err])
a: Shall be a rank-2, real or complex array containing the coefficient matrix.
It is an intent(in) argument.
pinva: Shall be a rank-2 array of the same kind as a, and size equal to that of transpose(a).
On output, it contains the Moore-Penrose pseudo-inverse of a.
rtol (optional): Shall be a scalar real value specifying the relative tolerance for singular value cutoff.
If not provided, the default threshold is , where is the
machine precision for the element type of a.
err (optional): Shall be a type(linalg_state_type) value. It is an intent(out) argument.
Computes the Moore-Penrose pseudo-inverse of the matrix , , and returns it in matrix pinva.
Raises LINALG_ERROR if the underlying SVD did not converge.
Raises LINALG_VALUE_ERROR if pinva and a have degenerate or incompatible sizes.
If err is not present, exceptions trigger an error stop.
! Matrix pseudo-inversion example: function, subroutine, and operator interfaces
program example_pseudoinverse
use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type
implicit none
real :: A(15,5), Am1(5,15)
type(linalg_state_type) :: state
! Generate random matrix A (15x15)
call random_number(A)
! Pseudo-inverse: Function interfcae
Am1 = pinv(A, err=state)
print *, 'Max error (function) : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! User threshold
Am1 = pinv(A, rtol=0.001, err=state)
print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! Pseudo-inverse: Subroutine interface
call pseudoinvert(A, Am1, err=state)
print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! Operator interface
Am1 = .pinv.A
print *, 'Max error (operator) : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
end program example_pseudoinverse
.pinv. - Moore-Penrose Pseudo-Inverse operatorExperimental
This operator returns the Moore-Penrose pseudo-inverse of a real or complex matrix .
The pseudo-inverse is computed using Singular Value Decomposition (SVD), and singular values
below a given threshold are treated as zero.
This interface is equivalent to the function pinv.
b = [[stdlib_linalg(module):operator(.pinv.)(interface)]] a
a: Shall be a rank-2 array of any real or complex kinds, with arbitrary dimensions . It is an intent(in) argument.
Returns a rank-2 array with the same type, kind, and rank as a, that contains the Moore-Penrose pseudo-inverse of a.
If an exception occurs, or if the input matrix is degenerate (e.g., rank-deficient), the returned matrix will contain NaNs.
For more detailed error handling, it is recommended to use the subroutine or function interfaces.
! Matrix pseudo-inversion example: function, subroutine, and operator interfaces
program example_pseudoinverse
use stdlib_linalg, only: pinv, pseudoinvert, operator(.pinv.), linalg_state_type
implicit none
real :: A(15,5), Am1(5,15)
type(linalg_state_type) :: state
! Generate random matrix A (15x15)
call random_number(A)
! Pseudo-inverse: Function interfcae
Am1 = pinv(A, err=state)
print *, 'Max error (function) : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! User threshold
Am1 = pinv(A, rtol=0.001, err=state)
print *, 'Max error (rtol=0.001): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! Pseudo-inverse: Subroutine interface
call pseudoinvert(A, Am1, err=state)
print *, 'Max error (subroutine): ', maxval(abs(A-matmul(A, matmul(Am1,A))))
! Operator interface
Am1 = .pinv.A
print *, 'Max error (operator) : ', maxval(abs(A-matmul(A, matmul(Am1,A))))
end program example_pseudoinverse
get_norm - Computes the vector norm of a generic-rank array.Experimental
This pure subroutine interface computes one of several vector norms of real or complex array , depending on
the value of the order input argument. may be an array of any rank.
Result nrm returns a real, scalar norm value for the whole array; if dim is specified, nrm is a rank n-1
array with the same shape as and dimension dim dropped, containing all norms evaluated along dim.
call get_norm (a, nrm, order, [, dim, err])
a: Shall be a rank-n real or complex array containing the data. It is an intent(in) argument.
nrm: if dim is absent, shall be a scalar with the norm evaluated over all the elements of the array. Otherwise, an array of rank n-1, and a shape similar
to that of a with dimension dim dropped.
order: Shall be an integer value or a character flag that specifies the norm type, as follows. It is an intent(in) argument.
| Integer input | Character Input | Norm type |
|---|---|---|
-huge(0) |
'-inf', '-Inf' |
Minimum absolute value ( \min_i{ \left |
1 |
'1' |
1-norm ( \sum_i{ \left |
2 |
'2' |
Euclidean norm |
>=3 |
'3','4',... |
p-norm ( \left( \sum_i{ \left |
huge(0) |
'inf', 'Inf' |
Maximum absolute value ( \max_i{ \left |
dim (optional): Shall be a scalar integer value with a value in the range from 1 to n, where n is the rank of the array. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument. If err is not present, the function is pure.
By default, the return value nrm is a scalar, and contains the norm as evaluated over all elements of the generic-rank array .
If the optional dim argument is present, nrm is a rank n-1 array with the same shape as except
for dimension dim, that is collapsed. Each element of nrm contains the 1D norm of the elements of ,
evaluated along dimension dim only.
Raises LINALG_ERROR if the requested norm type is invalid.
Raises LINALG_VALUE_ERROR if any of the arguments has an invalid size.
If err is not present, exceptions trigger an error stop.
! Vector norm: demonstrate usage of the function interface
program example_get_norm
use stdlib_linalg, only: get_norm, linalg_state_type
implicit none
real :: a(3,3), nrm, nrmd(3)
integer :: j
type(linalg_state_type) :: err
! a = [ -3.00000000 0.00000000 3.00000000
! -2.00000000 1.00000000 4.00000000
! -1.00000000 2.00000000 5.00000000 ]
a = reshape([(j-4,j=1,9)], [3,3])
print "(' a = [ ',3(g0,3x),2(/9x,3(g0,3x)),']')", transpose(a)
! Norm with integer input
call get_norm(a, nrm, 2)
print *, 'Euclidean norm = ',nrm ! 8.30662346
! Norm with character input
call get_norm(a, nrm, '2')
print *, 'Euclidean norm = ',nrm ! 8.30662346
! Euclidean norm of row arrays, a(i,:)
call get_norm(a, nrmd, 2, dim=2)
print *, 'Rows norms = ',nrmd ! 4.24264050 4.58257580 5.47722578
! Euclidean norm of columns arrays, a(:,i)
call get_norm(a, nrmd, 2, dim=1)
print *, 'Columns norms = ',nrmd ! 3.74165750 2.23606801 7.07106781
! Infinity norms
call get_norm(a, nrm, 'inf')
print *, 'maxval(||a||) = ',nrm ! 5.00000000
call get_norm(a, nrmd, 'inf', dim=2)
print *, 'maxval(||a(i,:)||) = ',nrmd ! 3.00000000 4.00000000 5.00000000
call get_norm(a, nrm, '-inf')
print *, 'minval(||a||) = ',nrm ! 0.00000000
call get_norm(a, nrmd, '-inf', dim=1)
print *, 'minval(||a(:,i)||) = ',nrmd ! 1.00000000 0.00000000 3.00000000
! Catch Error:
! [norm] returned Value Error: dimension 4 is out of rank for shape(a)= [3, 3]
call get_norm(a, nrmd, 'inf', dim=4, err=err)
print *, 'invalid: ',err%print()
end program example_get_norm
norm - Computes the vector norm of a generic-rank array.Experimental
This function computes one of several vector norms of real or complex array , depending on
the value of the order input argument. may be an array of any rank.
x = norm (a, order, [, dim, err])
a: Shall be a rank-n real or complex array containing the data. It is an intent(in) argument.
order: Shall be an integer value or a character flag that specifies the norm type, as follows. It is an intent(in) argument.
| Integer input | Character Input | Norm type |
|---|---|---|
-huge(0) |
'-inf', '-Inf' |
Minimum absolute value ( \min_i{ \left |
1 |
'1' |
1-norm ( \sum_i{ \left |
2 |
'2' |
Euclidean norm |
>=3 |
'3','4',... |
p-norm ( \left( \sum_i{ \left |
huge(0) |
'inf', 'Inf' |
Maximum absolute value ( \max_i{ \left |
dim (optional): Shall be a scalar integer value with a value in the range from 1 to n, where n is the rank of the array. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument. If err is not present, the function is pure.
By default, the return value x is a scalar, and contains the norm as evaluated over all elements of the generic-rank array .
If the optional dim argument is present, x is a rank n-1 array with the same shape as except
for dimension dim, that is dropped. Each element of x contains the 1D norm of the elements of ,
evaluated along dimension dim only.
Raises LINALG_ERROR if the requested norm type is invalid.
Raises LINALG_VALUE_ERROR if any of the arguments has an invalid size.
If err is not present, exceptions trigger an error stop.
! Vector norm: demonstrate usage of the function interface
program example_norm
use stdlib_linalg, only: norm, linalg_state_type
implicit none
real :: a(3,3),na
integer :: j
type(linalg_state_type) :: err
! a = [ -3.00000000 0.00000000 3.00000000
! -2.00000000 1.00000000 4.00000000
! -1.00000000 2.00000000 5.00000000 ]
a = reshape([(j-4,j=1,9)], [3,3])
print "(' a = [ ',3(g0,3x),2(/9x,3(g0,3x)),']')", transpose(a)
! Norm with integer input
print *, 'Euclidean norm = ',norm(a, 2) ! 8.30662346
! Norm with character input
print *, 'Euclidean norm = ',norm(a, '2') ! 8.30662346
! Euclidean norm of row arrays, a(i,:)
print *, 'Rows norms = ',norm(a, 2, dim=2) ! 4.24264050 4.58257580 5.47722578
! Euclidean norm of columns arrays, a(:,i)
print *, 'Columns norms = ',norm(a, 2, dim=1) ! 3.74165750 2.23606801 7.07106781
! Infinity norms
print *, 'maxval(||a||) = ',norm(a, 'inf') ! 5.00000000
print *, 'maxval(||a(i,:)||) = ',norm(a, 'inf', dim=2) ! 3.00000000 4.00000000 5.00000000
print *, 'minval(||a||) = ',norm(a, '-inf') ! 0.00000000
print *, 'minval(||a(:,i)||) = ',norm(a, '-inf', dim=1) ! 1.00000000 0.00000000 3.00000000
! Catch Error:
! [norm] returned Value Error: dimension 4 is out of rank for shape(a)= [3, 3]
print *, 'invalid: ',norm(a,'inf', dim=4, err=err)
print *, 'error = ',err%print()
end program example_norm
mnorm - Computes the matrix norm of a generic-rank array.Experimental
This function computes one of several matrix norms of real or complex array , depending on
the value of the order input argument. must be an array of rank 2 or higher. For arrays of rank > 2,
matrix norms are computed over specified dimensions.
x = mnorm (a [, order, dim, err])
a: Shall be a rank-n real or complex array containing the data, where n >= 2. It is an intent(in) argument.
order (optional): Shall be an integer value or a character flag that specifies the norm type, as follows. It is an intent(in) argument.
| Integer input | Character Input | Norm type |
|---|---|---|
1 |
'1' |
1-norm (maximum column sum) ( \max_j \sum_i{ \left |
2 |
'2' |
2-norm (largest singular value) |
| (not prov.) | 'Euclidean','Frobenius','Fro' |
Frobenius norm ( \sqrt{\sum_{i,j}{ \left |
huge(0) |
'inf', 'Inf', 'INF' |
Infinity norm (maximum row sum) ( \max_i \sum_j{ \left |
dim (optional): For arrays of rank > 2, shall be an integer array of size 2 specifying the dimensions over which to compute the matrix norm. Default value is [1,2]. It is an intent(in) argument.
err (optional): Shall be a type(linalg_state_type) value. This is an intent(out) argument.
For rank-2 input arrays, the return value x is a scalar containing the matrix norm.
For arrays of rank > 2, if the optional dim argument is present, x is a rank n-2 array with the same shape as except
for dimensions dim(1) and dim(2), which are dropped. Each element of x contains the matrix norm of the corresponding submatrix of ,
evaluated over the specified dimensions only, with the given order.
If an invalid norm type is provided, defaults to 1-norm and raises LINALG_ERROR.
Raises LINALG_VALUE_ERROR if any of the arguments has an invalid size.
If err is not present, exceptions trigger an error stop.
program example_mnorm
use stdlib_linalg, only: mnorm
use stdlib_kinds, only: sp
implicit none
real(sp) :: a(3,3), na
real(sp) :: b(3,3,4), nb(4) ! Array of 4 3x3 matrices
! Initialize example matrix
a = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
! Compute Euclidean norm of single matrix
na = mnorm(a, 'Euclidean')
print *, "Euclidean norm of matrix a:", na
! Initialize array of matrices
b(:,:,1) = a
b(:,:,2) = 2*a
b(:,:,3) = 3*a
b(:,:,4) = 4*a
! Compute infinity norm of each 3x3 matrix in b
nb = mnorm(b, 'inf', dim=[1,2])
! 18.0000000 36.0000000 54.0000000 72.0000000
print *, "Infinity norms of matrices in b:", nb
end program example_mnorm