Addition of diag, eye, and trace

src/CMakeLists.txt

@@ -3,6 +3,8 @@
 # Create a list of the files to be preprocessed
 set(fppFiles
     stdlib_experimental_io.fypp
+    stdlib_experimental_linalg.fypp
+    stdlib_experimental_linalg_diag.fypp
     stdlib_experimental_optval.fypp
     stdlib_experimental_stats.fypp
     stdlib_experimental_stats_mean.fypp

src/Makefile.manual

@@ -2,6 +2,8 @@ SRC = f18estop.f90 \
       stdlib_experimental_ascii.f90 \
       stdlib_experimental_error.f90 \
       stdlib_experimental_io.f90 \
+      stdlib_experimental_linalg.f90 \
+      stdlib_experimental_linalg_diag.f90 \
       stdlib_experimental_kinds.f90 \
       stdlib_experimental_optval.f90 \
       stdlib_experimental_quadrature.f90 \
@@ -42,6 +44,7 @@ stdlib_experimental_io.o: \
 	stdlib_experimental_error.o \
 	stdlib_experimental_optval.o \
 	stdlib_experimental_kinds.o
+stdlib_experimental_linalg_diag.o: stdlib_experimental_kinds.o
 stdlib_experimental_optval.o: stdlib_experimental_kinds.o
 stdlib_experimental_quadrature.o: stdlib_experimental_kinds.o
 stdlib_experimental_stats_mean.o: \
@@ -59,6 +62,8 @@ stdlib_experimental_stats_var.o: \
 
 # Fortran sources that are built from fypp templates
 stdlib_experimental_io.f90: stdlib_experimental_io.fypp
+stdlib_experimental_linalg.f90: stdlib_experimental_linalg.fypp
+stdlib_experimental_linalg_diag.f90: stdlib_experimental_linalg_diag.fypp
 stdlib_experimental_quadrature.f90: stdlib_experimental_quadrature.fypp
 stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp
 stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp

src/stdlib_experimental_linalg.fypp

@@ -0,0 +1,80 @@
+#:include "common.fypp"
+#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
+module stdlib_experimental_linalg
+  use stdlib_experimental_kinds, only: sp, dp, qp, &
+    int8, int16, int32, int64
+  implicit none
+  private
+
+  public :: diag
+  public :: eye
+  public :: trace
+
+  interface diag
+      !
+      ! Vector to matrix
+      !
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module function diag_${t1[0]}$${k1}$(v) result(res)
+        ${t1}$, intent(in) :: v(:)
+        ${t1}$ :: res(size(v),size(v))
+      end function diag_${t1[0]}$${k1}$
+    #:endfor
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module function diag_${t1[0]}$${k1}$_k(v,k) result(res)
+        ${t1}$, intent(in) :: v(:)
+        integer, intent(in) :: k
+        ${t1}$ :: res(size(v)+abs(k),size(v)+abs(k))
+      end function diag_${t1[0]}$${k1}$_k
+    #:endfor
+
+      !
+      ! Matrix to vector
+      !
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module function diag_${t1[0]}$${k1}$_mat(A) result(res)
+        ${t1}$, intent(in) :: A(:,:)
+        ${t1}$ :: res(minval(shape(A)))
+      end function diag_${t1[0]}$${k1}$_mat
+    #:endfor
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
+        ${t1}$, intent(in) :: A(:,:)
+        integer, intent(in) :: k
+        ${t1}$ :: res(minval(shape(A))-abs(k))
+      end function diag_${t1[0]}$${k1}$_mat_k
+    #:endfor
+  end interface
+
+  ! Matrix trace
+  interface trace
+    #:for k1, t1 in RCI_KINDS_TYPES
+      module procedure trace_${t1[0]}$${k1}$
+    #:endfor
+  end interface
+
+contains
+
+    function eye(n) result(res)
+      integer, intent(in) :: n
+      integer(int8) :: res(n, n)
+      integer :: i
+      res = 0
+      do i = 1, n
+         res(i, i) = 1
+      end do
+    end function eye
+
+
+    #:for k1, t1 in RCI_KINDS_TYPES
+      function trace_${t1[0]}$${k1}$(A) result(res)
+        ${t1}$, intent(in) :: A(:,:)
+        ${t1}$ :: res
+        integer :: i
+        res = 0
+        do i = 1, minval(shape(A))
+          res = res + A(i,i)
+        end do
+      end function trace_${t1[0]}$${k1}$
+    #:endfor
+end module

src/stdlib_experimental_linalg.md

@@ -0,0 +1,156 @@
+# Linear Algebra
+
+* [`diag` - Create a diagonal array or extract the diagonal elements of an array](#diag---create-a-diagonal-array-or-extract-the-diagonal-elements-of-an-array)
+* [`eye` - Construct the identity matrix](#eye---construct-the-identity-matrix)
+* [`trace` - Trace of a matrix](#trace---trace-of-a-matrix)
+
+## `diag` - Create a diagonal array or extract the diagonal elements of an array
+
+### Description
+
+Create a diagonal array or extract the diagonal elements of an array
+
+### Syntax
+
+`d = diag(a [, k])`
+
+### Arguments
+
+`a`: Shall be a rank-1 or or rank-2 array. If `a` is a rank-1 array (i.e. a vector) then `diag` returns a rank-2 array with the elements of `a` on the diagonal. If `a` is a rank-2 array (i.e. a matrix) then `diag` returns a rank-1 array of the diagonal elements.
+
+`k` (optional): Shall be a scalar of type `integer` and specifies the diagonal. The default `k = 0` represents the main diagonal, `k > 0` are diagonals above the main diagonal, `k < 0` are diagonals below the main diagonal.
+
+### Return value
+
+Returns a diagonal array or a vector with the extracted diagonal elements.
+
+### Example
+
+```fortran
+program demo_diag1
+    use stdlib_experimental_linalg, only: diag
+    implicit none
+    real, allocatable :: A(:,:)
+    integer :: i
+    A = diag([(1,i=1,10)]) ! creates a 10 by 10 identity matrix
+end program demo_diag1
+```
+
+```fortran
+program demo_diag2
+    use stdlib_experimental_linalg, only: diag
+    implicit none
+    real :: v(:)
+    real, allocatable :: A(:,:)
+    integer :: i
+    v = [1,2,3,4,5]
+    A = diag(v) ! creates a 5 by 5 matrix with elements of v on the diagonal
+end program demo_diag2
+```
+
+```fortran
+program demo_diag3
+    use stdlib_experimental_linalg, only: diag
+    implicit none
+    integer, parameter :: n = 10
+    real :: c(n), ul(n-1)
+    real :: A(n,n)
+    integer :: i
+    c = 2
+    ul = -1
+    A = diag(ul,-1) + diag(c) + diag(ul,1) ! Gil Strang's favorite matrix
+end program demo_diag3
+```
+
+```fortran
+program demo_diag4
+    use stdlib_experimental_linalg, only: diag
+    implicit none
+    integer, parameter :: n = 12
+    real :: A(n,n)
+    real :: v(n)
+    integer :: i
+    call random_number(A)
+    v = diag(A) ! v contains diagonal elements of A
+end program demo_diag4
+```
+
+```fortran
+program demo_diag5
+    use stdlib_experimental_linalg, only: diag
+    implicit none
+    integer, parameter :: n = 3
+    real :: A(n,n)
+    real, allocatable :: v(:)
+    integer :: i
+    A = reshape([1,2,3,4,5,6,7,8,9],[n,n])
+    v = diag(A,-1) ! v is [2,6]
+    v = diag(A,1)  ! v is [4,8]
+end program demo_diag5
+```
+
+## `eye` - Construct the identity matrix
+
+### Description
+
+Construct the identity matrix
+
+## Syntax
+
+`I = eye(n)`
+
+### Arguments
+
+`n`: Shall be a scalar of default type `integer`. 
+
+### Return value
+
+Returns the identity matrix, i.e. a square matrix with ones on the main diagonal and zeros elsewhere. The return value is of type `integer(int8)`.
+
+### Example
+
+```fortran
+program demo_eye1
+    use stdlib_experimental_linalg, only: eye
+    implicit none
+    real :: a(3,3)
+    A = eye(3)
+end program demo_eye1
+```
+
+```fortran
+program demo_eye2
+    use stdlib_experimental_linalg, only: eye, diag
+    implicit none
+    print *, all(eye(4) == diag([1,1,1,1])) ! prints .true.
+end program demo_eye2
+```
+
+## `trace` - Trace of a matrix
+
+### Description
+
+Trace of a matrix (rank-2 array)
+
+### Syntax
+
+`result = trace(A)`
+
+### Arguments
+
+`A`: Shall be a rank-2 array. If `A` is not square, then `trace(A)` will return the sum of diagonal values from the square sub-section of `A`.
+
+### Return value
+
+Returns the trace of the matrix, i.e. the sum of diagonal elements.
+
+### Example
+```fortran
+program demo_trace
+    use stdlib_experimental_linalg, only: trace
+    implicit none
+    real :: A(3,3)
+    A = reshape([1,2,3,4,5,6,7,8,9],[3,3])
+    print *, trace(A) ! 1 + 5 + 9
+end program demo_trace
+```

src/stdlib_experimental_linalg_diag.fypp

@@ -0,0 +1,80 @@
+#:include "common.fypp"
+#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
+submodule (stdlib_experimental_linalg) stdlib_experimental_linalg_diag
+
+  implicit none
+
+contains
+
+    #:for k1, t1 in RCI_KINDS_TYPES
+      function diag_${t1[0]}$${k1}$(v) result(res)
+        ${t1}$, intent(in) :: v(:)
+        ${t1}$ :: res(size(v),size(v))
+        integer :: i
+        res = 0
+        do i = 1, size(v)
+          res(i,i) = v(i)
+        end do
+      end function diag_${t1[0]}$${k1}$
+    #:endfor
+
+
+    #:for k1, t1 in RCI_KINDS_TYPES
+      function diag_${t1[0]}$${k1}$_k(v,k) result(res)
+        ${t1}$, intent(in) :: v(:)
+        integer, intent(in) :: k
+        ${t1}$ :: res(size(v)+abs(k),size(v)+abs(k))
+        integer :: i, sz
+        sz = size(v)
+        res = 0
+        if (k > 0) then
+          do i = 1, sz
+              res(i,k+i) = v(i)
+          end do
+        else if (k < 0) then
+          do i = 1, sz
+              res(i+abs(k),i) = v(i)
+          end do
+        else
+          do i = 1, sz
+              res(i,i) = v(i)
+          end do
+        end if
+      end function diag_${t1[0]}$${k1}$_k
+    #:endfor
+
+    #:for k1, t1 in RCI_KINDS_TYPES
+      function diag_${t1[0]}$${k1}$_mat(A) result(res)
+        ${t1}$, intent(in) :: A(:,:)
+        ${t1}$ :: res(minval(shape(A)))
+        integer :: i
+        do i = 1, minval(shape(A))
+          res(i) = A(i,i)
+        end do
+      end function diag_${t1[0]}$${k1}$_mat
+    #:endfor
+
+    #:for k1, t1 in RCI_KINDS_TYPES
+      function diag_${t1[0]}$${k1}$_mat_k(A,k) result(res)
+        ${t1}$, intent(in) :: A(:,:)
+        integer, intent(in) :: k
+        ${t1}$ :: res(minval(shape(A))-abs(k))
+        integer :: i, sz
+        sz = minval(shape(A))-abs(k)
+        if (k > 0) then
+          do i = 1, sz
+              res(i) = A(i,k+i)
+          end do
+        else if (k < 0) then
+          do i = 1, sz
+              res(i) = A(i+abs(k),i)
+          end do
+        else
+          do i = 1, sz
+              res(i) = A(i,i)
+          end do
+        end if
+      end function diag_${t1[0]}$${k1}$_mat_k
+    #:endfor
+
+end submodule

src/tests/CMakeLists.txt

@@ -8,6 +8,7 @@ endmacro(ADDTEST)
 
 add_subdirectory(ascii)
 add_subdirectory(io)
+add_subdirectory(linalg)
 add_subdirectory(optval)
 add_subdirectory(stats)
 add_subdirectory(system)

src/tests/linalg/CMakeLists.txt

@@ -0,0 +1,2 @@
+ADDTEST(linalg)
+