Merge branch 'fortran-lang:master' into subprocess

Author: perazz

doc/specs/stdlib_linalg.md

@@ -239,37 +239,34 @@ Pure function.
 
 ### Description
 
-Construct the identity matrix.
+Constructs the identity matrix.
 
 ### Syntax
 
-`I = ` [[stdlib_linalg(module):eye(function)]] `(dim1 [, dim2])`
+`I = ` [[stdlib_linalg(module):eye(function)]] `(dim1 [, dim2] [, mold])`
 
 ### Arguments
 
-`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. 
+- `dim1`: A scalar of type `integer`. This is an `intent(in)` argument and specifies the number of rows.
+- `dim2`: A scalar of type `integer`. This is an optional `intent(in)` argument specifying the number of columns. If not provided, the matrix is square (`dim1 = dim2`).
+- `mold`: A scalar of any supported `integer`, `real`, or `complex` type. This is an optional `intent(in)` argument. If provided, the returned identity matrix will have the same type and kind as `mold`. If not provided, the matrix will be of type `real(real64)` by default.
 
 ### Return value
 
-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.
+Returns the identity matrix, with ones on the main diagonal and zeros elsewhere. 
+
+- By default, the return value is of type `real(real64)`, which is recommended for arithmetic safety.
+- If the `mold` argument is provided, the return value will match the type and kind of `mold`, allowing for arbitrary `integer`, `real`, or `complex` return types.
 
-#### Warning
+### Example
 
-Since the result of `eye` is of `integer(int8)` type, one should be careful about using it in arithmetic expressions. For example:
 ```fortran
-!> Be careful
-A = eye(2,2)/2     !! A == 0.0
-!> Recommend
-A = eye(2,2)/2.0   !! A == diag([0.5, 0.5])
+!> Return default type (real64)
+A = eye(2,2)/2             !! A == diag([0.5_dp, 0.5_dp])
+!> Return 32-bit complex
+A = eye(2,2, mold=(0.0,0.0))/2   !! A == diag([(0.5,0.5), (0.5,0.5)])
 ```
 
-### Example
-
 ```fortran
 {!example/linalg/example_eye1.f90!}
 ```
@@ -509,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_eye1.f90

@@ -5,10 +5,10 @@ program example_eye1
   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)]
+  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

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
@@ -189,6 +190,16 @@ module stdlib_linalg
     #:endfor
   end interface
 
+  ! Identity matrix 
+  interface eye
+    !! version: experimental
+    !!
+    !! Constructs the identity matrix
+    !! ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))    
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module procedure eye_${t1[0]}$${k1}$
+    #:endfor
+  end interface eye
 
   ! Outer product (of two vectors)
   interface outer_product
@@ -299,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
@@ -1411,24 +1447,27 @@ contains
     !>
     !> Constructs the identity matrix.
     !> ([Specification](../page/specs/stdlib_linalg.html#eye-construct-the-identity-matrix))
-    pure function eye(dim1, dim2) result(result)
+    #:for k1, t1 in RCI_KINDS_TYPES
+    pure function eye_${t1[0]}$${k1}$(dim1, dim2, mold) result(result)
 
         integer, intent(in) :: dim1
         integer, intent(in), optional :: dim2
-        integer(int8), allocatable :: result(:, :)
+        ${t1}$, intent(in) #{if t1 == 'real(dp)'}#, optional #{endif}#:: mold        
+        ${t1}$, allocatable :: result(:, :)
 
         integer :: dim2_
         integer :: i
 
         dim2_ = optval(dim2, dim1)
         allocate(result(dim1, dim2_))
 
-        result = 0_int8
+        result = 0
         do i = 1, min(dim1, dim2_)
-            result(i, i) = 1_int8
+            result(i, i) = 1
         end do
 
-    end function eye
+    end function eye_${t1[0]}$${k1}$
+    #:endfor
 
     #:for k1, t1 in RCI_KINDS_TYPES
       function trace_${t1[0]}$${k1}$(A) result(res)
@@ -1555,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