diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 2f03f0fad..f1ed82b02 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -509,6 +509,37 @@ Returns a `logical` scalar that is `.true.` if the input matrix is skew-symmetri {!example/linalg/example_is_skew_symmetric.f90!} ``` +## `hermitian` - Compute the Hermitian version of a rank-2 matrix + +### Status + +Experimental + +### Description + +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)`). + +### Syntax + +`h = ` [[stdlib_linalg(module):hermitian(interface)]] `(a)` + +### Arguments + +`a`: Shall be a rank-2 array of type `integer`, `real`, or `complex`. The input matrix `a` is not modified. + +### Return value + +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 `transpose(a)`. + +### Example + +```fortran +{!example/linalg/example_hermitian.f90!} +``` + ## `is_hermitian` - Checks if a matrix is Hermitian ### Status diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index fc3f823d9..27f5ccbc5 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -6,6 +6,7 @@ ADD_EXAMPLE(diag5) ADD_EXAMPLE(eye1) ADD_EXAMPLE(eye2) ADD_EXAMPLE(is_diagonal) +ADD_EXAMPLE(hermitian) ADD_EXAMPLE(is_hermitian) ADD_EXAMPLE(is_hessenberg) ADD_EXAMPLE(is_skew_symmetric) diff --git a/example/linalg/example_hermitian.f90 b/example/linalg/example_hermitian.f90 new file mode 100644 index 000000000..234e4b7b6 --- /dev/null +++ b/example/linalg/example_hermitian.f90 @@ -0,0 +1,19 @@ +! 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 diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index ad36ad677..663a0d5ae 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -48,6 +48,7 @@ module stdlib_linalg public :: is_diagonal public :: is_symmetric public :: is_skew_symmetric + public :: hermitian public :: is_hermitian public :: is_triangular public :: is_hessenberg @@ -298,6 +299,31 @@ module stdlib_linalg #:endfor end interface is_hermitian + interface hermitian + !! version: experimental + !! + !! Computes the Hermitian version of a rank-2 matrix. + !! For complex matrices, this returns `conjg(transpose(a))`. + !! For real or integer matrices, this returns `transpose(a)`. + !! + !! Usage: + !! ``` + !! A = reshape([(1, 2), (3, 4), (5, 6), (7, 8)], [2, 2]) + !! AH = hermitian(A) + !! ``` + !! + !! [Specification](../page/specs/stdlib_linalg.html#hermitian-compute-the-hermitian-version-of-a-rank-2-matrix) + !! + + #:for k1, t1 in RCI_KINDS_TYPES + pure module function hermitian_${t1[0]}$${k1}$(a) result(ah) + ${t1}$, intent(in) :: a(:,:) + ${t1}$ :: ah(size(a, 2), size(a, 1)) + end function hermitian_${t1[0]}$${k1}$ + #:endfor + + end interface hermitian + ! Check for triangularity interface is_triangular @@ -1471,6 +1497,17 @@ contains end function is_hermitian_${t1[0]}$${k1}$ #:endfor + #:for k1, t1 in RCI_KINDS_TYPES + pure module function hermitian_${t1[0]}$${k1}$(a) result(ah) + ${t1}$, intent(in) :: a(:,:) + ${t1}$ :: ah(size(a, 2), size(a, 1)) + #:if t1.startswith('complex') + ah = conjg(transpose(a)) + #:else + ah = transpose(a) + #:endif + end function hermitian_${t1[0]}$${k1}$ + #:endfor #:for k1, t1 in RCI_KINDS_TYPES function is_triangular_${t1[0]}$${k1}$(A,uplo) result(res) @@ -1546,4 +1583,5 @@ contains end function is_hessenberg_${t1[0]}$${k1}$ #:endfor + end module stdlib_linalg diff --git a/test/linalg/test_linalg.fypp b/test/linalg/test_linalg.fypp index fcb9c6abb..ef4e89cf7 100644 --- a/test/linalg/test_linalg.fypp +++ b/test/linalg/test_linalg.fypp @@ -4,7 +4,7 @@ module test_linalg use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64 - use stdlib_linalg, only: diag, eye, trace, outer_product, cross_product, kronecker_product + use stdlib_linalg, only: diag, eye, trace, outer_product, cross_product, kronecker_product, hermitian use stdlib_linalg_state, only: linalg_state_type, LINALG_SUCCESS, linalg_error_handling implicit none @@ -52,6 +52,9 @@ contains new_unittest("trace_int64", test_trace_int64), & #:for k1, t1 in RCI_KINDS_TYPES new_unittest("kronecker_product_${t1[0]}$${k1}$", test_kronecker_product_${t1[0]}$${k1}$), & + #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + new_unittest("hermitian_${t1[0]}$${k1}$", test_hermitian_${t1[0]}$${k1}$), & #:endfor new_unittest("outer_product_rsp", test_outer_product_rsp), & new_unittest("outer_product_rdp", test_outer_product_rdp), & @@ -597,6 +600,32 @@ contains end subroutine test_kronecker_product_${t1[0]}$${k1}$ #:endfor + #:for k1, t1 in CMPLX_KINDS_TYPES + subroutine test_hermitian_${t1[0]}$${k1}$(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + integer, parameter :: m = 2, n = 3 + ${t1}$, dimension(m,n) :: A + ${t1}$, dimension(n,m) :: AT, expected, diff + real(${k1}$), parameter :: tol = 1.e-6_${k1}$ + + integer :: i,j + + do concurrent (i=1:m,j=1:n) + A (i,j) = cmplx(i,-j,kind=${k1}$) + expected(j,i) = cmplx(i,+j,kind=${k1}$) + end do + + + AT = hermitian(A) + + diff = AT - expected + + call check(error, all(abs(diff) < abs(tol)), "hermitian: all(abs(diff) < abs(tol)) failed") + + end subroutine test_hermitian_${t1[0]}$${k1}$ + #:endfor + subroutine test_outer_product_rsp(error) !> Error handling type(error_type), allocatable, intent(out) :: error