diag - Create a diagonal array or extract the diagonal elements of an arrayExperimental
Create a diagonal array or extract the diagonal elements of an array
d = diag(a [, k])
a: Shall be a rank-1 or or rank-2 array. If a is a rank-1 array (i.e. a vector) then diag returns a rank-2 array with the elements of a on the diagonal. If a is a rank-2 array (i.e. a matrix) then diag returns a rank-1 array of the diagonal elements.
k (optional): Shall be a scalar of type integer and specifies the diagonal. The default k = 0 represents the main diagonal, k > 0 are diagonals above the main diagonal, k < 0 are diagonals below the main diagonal.
Returns a diagonal array or a vector with the extracted diagonal elements.
program example_diag1
use stdlib_linalg, only: diag
implicit none
real, allocatable :: A(:, :)
integer :: i
A = diag([(1, i=1, 10)]) ! creates a 10 by 10 identity matrix
end program example_diag1
program example_diag2
use stdlib_linalg, only: diag
implicit none
real, allocatable :: v(:)
real, allocatable :: A(:, :)
v = [1, 2, 3, 4, 5]
A = diag(v) ! creates a 5 by 5 matrix with elements of v on the diagonal
end program example_diag2
program example_diag3
use stdlib_linalg, only: diag
implicit none
integer, parameter :: n = 10
real :: c(n), ul(n - 1)
real :: A(n, n)
c = 2
ul = -1
A = diag(ul, -1) + diag(c) + diag(ul, 1) ! Gil Strang's favorite matrix
end program example_diag3
program example_diag4
use stdlib_linalg, only: diag
implicit none
integer, parameter :: n = 12
real :: A(n, n)
real :: v(n)
call random_number(A)
v = diag(A) ! v contains diagonal elements of A
end program example_diag4
program example_diag5
use stdlib_linalg, only: diag
implicit none
integer, parameter :: n = 3
real :: A(n, n)
real, allocatable :: v(:)
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [n, n])
v = diag(A, -1) ! v is [2,6]
v = diag(A, 1) ! v is [4,8]
end program example_diag5
eye - Construct the identity matrixExperimental
Pure function.
Construct the identity matrix.
I = eye(dim1 [, dim2])
dim1: Shall be a scalar of default type integer.
This is an intent(in) argument.
dim2: Shall be a scalar of default type integer.
This is an intent(in) and optional argument.
Return the identity matrix, i.e. a matrix with ones on the main diagonal and zeros elsewhere. The return value is of type integer(int8).
The use of int8 was suggested to save storage.
Since the result of eye is of integer(int8) type, one should be careful about using it in arithmetic expressions. For example:
!> Be careful
A = eye(2,2)/2 !! A == 0.0
!> Recommend
A = eye(2,2)/2.0 !! A == diag([0.5, 0.5])
program example_eye1
use stdlib_linalg, only: eye
implicit none
integer :: i(2, 2)
real :: a(3, 3)
real :: b(2, 3) !! Matrix is non-square.
complex :: c(2, 2)
I = eye(2) !! [1,0; 0,1]
A = eye(3) !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
A = eye(3, 3) !! [1.0,0.0,0.0; 0.0,1.0,0.0; 0.0,0.0,1.0]
B = eye(2, 3) !! [1.0,0.0,0.0; 0.0,1.0,0.0]
C = eye(2, 2) !! [(1.0,0.0),(0.0,0.0); (0.0,0.0),(1.0,0.0)]
C = (1.0, 1.0)*eye(2, 2) !! [(1.0,1.0),(0.0,0.0); (0.0,0.0),(1.0,1.0)]
end program example_eye1
program example_eye2
use stdlib_linalg, only: eye, diag
implicit none
print *, all(eye(4) == diag([1, 1, 1, 1])) ! prints .true.
end program example_eye2
trace - Trace of a matrixExperimental
Trace of a matrix (rank-2 array)
result = trace(A)
A: Shall be a rank-2 array. If A is not square, then trace(A) will return the sum of diagonal values from the square sub-section of A.
Returns the trace of the matrix, i.e. the sum of diagonal elements.
program example_trace
use stdlib_linalg, only: trace
implicit none
real :: A(3, 3)
A = reshape([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 3])
print *, trace(A) ! 1 + 5 + 9
end program example_trace
outer_product - Computes the outer product of two vectorsExperimental
Computes the outer product of two vectors
d = outer_product(u, v)
u: Shall be a rank-1 array
v: Shall be a rank-1 array
Returns a rank-2 array equal to u v^T (where u, v are considered column vectors). The shape of the returned array is [size(u), size(v)].
program example_outer_product
use stdlib_linalg, only: outer_product
implicit none
real, allocatable :: A(:, :), u(:), v(:)
u = [1., 2., 3.]
v = [3., 4.]
A = outer_product(u, v)
!A = reshape([3., 6., 9., 4., 8., 12.], [3,2])
end program example_outer_product
kronecker_product - Computes the Kronecker product of two rank-2 arraysExperimental
Computes the Kronecker product of two rank-2 arrays
C = kronecker_product(A, B)
A: Shall be a rank-2 array with dimensions M1, N1
B: Shall be a rank-2 array with dimensions M2, N2
Returns a rank-2 array equal to A \otimes B. The shape of the returned array is [M1*M2, N1*N2].
program example_kronecker_product
use stdlib_linalg, only: kronecker_product
implicit none
integer, parameter :: m1 = 1, n1 = 2, m2 = 2, n2 = 3
integer :: i, j
real :: A(m1, n1), B(m2,n2)
real, allocatable :: C(:,:)
do j = 1, n1
do i = 1, m1
A(i,j) = i*j ! A = [1, 2]
end do
end do
do j = 1, n2
do i = 1, m2 ! B = [ 1, 2, 3 ]
B(i,j) = i*j ! [ 2, 4, 6 ]
end do
end do
C = kronecker_product(A, B)
! C = [ a(1,1) * B(:,:) | a(1,2) * B(:,:) ]
! or in other words,
! C = [ 1.00 2.00 3.00 2.00 4.00 6.00 ]
! [ 2.00 4.00 6.00 4.00 8.00 12.00 ]
end program example_kronecker_product
cross_product - Computes the cross product of two vectorsExperimental
Computes the cross product of two vectors
c = cross_product(a, b)
a: Shall be a rank-1 and size-3 array
b: Shall be a rank-1 and size-3 array
Returns a rank-1 and size-3 array which is perpendicular to both a and b.
program demo_cross_product
use stdlib_linalg, only: cross_product
implicit none
real :: a(3), b(3), c(3)
a = [1., 0., 0.]
b = [0., 1., 0.]
c = cross_product(a, b)
!c = [0., 0., 1.]
end program demo_cross_product
is_square - Checks if a matrix is squareExperimental
Checks if a matrix is square
d = is_square(A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is square, and .false. otherwise.
program example_is_square
use stdlib_linalg, only: is_square
implicit none
real :: A(2, 2), B(3, 2)
logical :: res
A = reshape([1., 2., 3., 4.], shape(A))
B = reshape([1., 2., 3., 4., 5., 6.], shape(B))
res = is_square(A) ! returns .true.
res = is_square(B) ! returns .false.
end program example_is_square
is_diagonal - Checks if a matrix is diagonalExperimental
Checks if a matrix is diagonal
d = is_diagonal(A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is diagonal, and .false. otherwise.
Note that nonsquare matrices may be diagonal, so long as a_ij = 0 when i /= j.
program example_is_diagonal
use stdlib_linalg, only: is_diagonal
implicit none
real :: A(2, 2), B(2, 2)
logical :: res
A = reshape([1., 0., 0., 4.], shape(A))
B = reshape([1., 0., 3., 4.], shape(B))
res = is_diagonal(A) ! returns .true.
res = is_diagonal(B) ! returns .false.
end program example_is_diagonal
is_symmetric - Checks if a matrix is symmetricExperimental
Checks if a matrix is symmetric
d = is_symmetric(A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is symmetric, and .false. otherwise.
program example_is_symmetric
use stdlib_linalg, only: is_symmetric
implicit none
real :: A(2, 2), B(2, 2)
logical :: res
A = reshape([1., 3., 3., 4.], shape(A))
B = reshape([1., 0., 3., 4.], shape(B))
res = is_symmetric(A) ! returns .true.
res = is_symmetric(B) ! returns .false.
end program example_is_symmetric
is_skew_symmetric - Checks if a matrix is skew-symmetricExperimental
Checks if a matrix is skew-symmetric
d = is_skew_symmetric(A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is skew-symmetric, and .false. otherwise.
program example_is_skew_symmetric
use stdlib_linalg, only: is_skew_symmetric
implicit none
real :: A(2, 2), B(2, 2)
logical :: res
A = reshape([0., -3., 3., 0.], shape(A))
B = reshape([0., 3., 3., 0.], shape(B))
res = is_skew_symmetric(A) ! returns .true.
res = is_skew_symmetric(B) ! returns .false.
end program example_is_skew_symmetric
is_hermitian - Checks if a matrix is HermitianExperimental
Checks if a matrix is Hermitian
d = is_hermitian(A)
A: Shall be a rank-2 array
Returns a logical scalar that is .true. if the input matrix is Hermitian, and .false. otherwise.
program example_is_hermitian
use stdlib_linalg, only: is_hermitian
implicit none
complex :: A(2, 2), B(2, 2)
logical :: res
A = reshape([cmplx(1., 0.), cmplx(3., -1.), cmplx(3., 1.), cmplx(4., 0.)], shape(A))
B = reshape([cmplx(1., 0.), cmplx(3., 1.), cmplx(3., 1.), cmplx(4., 0.)], shape(B))
res = is_hermitian(A) ! returns .true.
res = is_hermitian(B) ! returns .false.
end program example_is_hermitian
is_triangular - Checks if a matrix is triangularExperimental
Checks if a matrix is triangular
d = is_triangular(A,uplo)
A: Shall be a rank-2 array
uplo: Shall be a single character from {'u','U','l','L'}
Returns a logical scalar that is .true. if the input matrix is the type of triangular specified by uplo (upper or lower), and .false. otherwise.
Note that the definition of triangular used in this implementation allows nonsquare matrices to be triangular.
Specifically, upper triangular matrices satisfy a_ij = 0 when j < i, and lower triangular matrices satisfy a_ij = 0 when j > i.
program example_is_triangular
use stdlib_linalg, only: is_triangular
implicit none
real :: A(3, 3), B(3, 3)
logical :: res
A = reshape([1., 0., 0., 4., 5., 0., 7., 8., 9.], shape(A))
B = reshape([1., 0., 3., 4., 5., 0., 7., 8., 9.], shape(B))
res = is_triangular(A, 'u') ! returns .true.
res = is_triangular(B, 'u') ! returns .false.
end program example_is_triangular
is_hessenberg - Checks if a matrix is hessenbergExperimental
Checks if a matrix is Hessenberg
d = is_hessenberg(A,uplo)
A: Shall be a rank-2 array
uplo: Shall be a single character from {'u','U','l','L'}
Returns a logical scalar that is .true. if the input matrix is the type of Hessenberg specified by uplo (upper or lower), and .false. otherwise.
Note that the definition of Hessenberg used in this implementation allows nonsquare matrices to be Hessenberg.
Specifically, upper Hessenberg matrices satisfy a_ij = 0 when j < i-1, and lower Hessenberg matrices satisfy a_ij = 0 when j > i+1.
program example_is_hessenberg
use stdlib_linalg, only: is_hessenberg
implicit none
real :: A(3, 3), B(3, 3)
logical :: res
A = reshape([1., 2., 0., 4., 5., 6., 7., 8., 9.], shape(A))
B = reshape([1., 2., 3., 4., 5., 6., 7., 8., 9.], shape(B))
res = is_hessenberg(A, 'u') ! returns .true.
res = is_hessenberg(B, 'u') ! returns .false.
end program example_is_hessenberg