Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

linalg: hermitian #896

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 31 additions & 0 deletions doc/specs/stdlib_linalg.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions example/linalg/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
19 changes: 19 additions & 0 deletions example/linalg/example_hermitian.f90
Original file line number Diff line number Diff line change
@@ -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
38 changes: 38 additions & 0 deletions src/stdlib_linalg.fypp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -1546,4 +1583,5 @@ contains
end function is_hessenberg_${t1[0]}$${k1}$
#:endfor


end module stdlib_linalg
31 changes: 30 additions & 1 deletion test/linalg/test_linalg.fypp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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), &
Expand Down Expand Up @@ -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
Expand Down
Loading