Skip to content

Commit

Permalink
linalg: hermitian (#896)
Browse files Browse the repository at this point in the history
  • Loading branch information
perazz authored Dec 26, 2024
2 parents fcb5a50 + 5efba9b commit 095219e
Show file tree
Hide file tree
Showing 5 changed files with 118 additions and 1 deletion.
31 changes: 31 additions & 0 deletions doc/specs/stdlib_linalg.md
Original file line number Diff line number Diff line change
Expand Up @@ -506,6 +506,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 the intrinsic `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
37 changes: 37 additions & 0 deletions src/stdlib_linalg.fypp
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,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 @@ -309,6 +310,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 @@ -1568,6 +1594,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
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

0 comments on commit 095219e

Please sign in to comment.