linalg: hermitian (#896)

Author: perazz

doc/specs/stdlib_linalg.md

@@ -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

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)

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

src/stdlib_linalg.fypp

@@ -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
@@ -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
@@ -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)

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