Added tests for diag

Author: ivan-pi

src/tests/linalg/test_linalg.f90

@@ -3,22 +3,8 @@ program test_linalg
   use stdlib_experimental_kinds, only: sp, dp, qp, int8, int16, int32, int64
   use stdlib_experimental_linalg, only: diag, eye, trace
   implicit none
-
-  real(sp) :: a(5)
-  real(sp), allocatable :: d(:,:)
-  integer :: i, j
   logical :: warn
-
-  warn = .true.
-  a = [1,2,3,4,5]
-  d = diag(a,0)
-
-  do i = 1, size(d,1)
-    write(*,*) (d(i,j),j=1,size(d,2))
-  end do
-
-  a = diag(d,0)
-  print *, "a = ", a
+  warn = .false.
 
   !
   ! eye
@@ -28,7 +14,19 @@ program test_linalg
   !
   ! diag
   !
+  call test_diag_rsp
+  call test_diag_rsp_k
+  call test_diag_rdp
+  call test_diag_rqp
+
+  call test_diag_csp
+  call test_diag_cdp
+  call test_diag_cqp
 
+  call test_diag_int8
+  call test_diag_int16
+  call test_diag_int32
+  call test_diag_int64
 
   !
   ! trace
@@ -69,9 +67,198 @@ subroutine test_eye
           msg="abs(trace(cye) - complex(7.0_sp,0.0_sp)) < epsilon(1.0_sp) failed.",warn=warn)
   end subroutine
 
+  subroutine test_diag_rsp
+    integer, parameter :: n = 3
+    real(sp) :: v(n), a(n,n), b(n,n)
+    integer :: i,j
+    write(*,*) "test_diag_rsp"
+    v = [(i,i=1,n)]
+    a = diag(v)
+    b = reshape([((merge(i,0,i==j), i=1,n), j=1,n)], [n,n])
+    call check(all(a == b), &
+            msg="all(a == b) failed.",warn=warn)
+
+    call check(all(diag(3*a) == 3*v), &
+            msg="all(diag(3*a) == 3*v) failed.",warn=warn)
+  end subroutine
+
+  subroutine test_diag_rsp_k
+    integer, parameter :: n = 4
+    real(sp) :: a(n,n), b(n,n)
+    integer :: i,j
+    write(*,*) "test_diag_rsp_k"
+
+    a = diag([(1._sp,i=1,n-1)],-1)
+    
+    b = reshape([((merge(1,0,i==j+1), i=1,n), j=1,n)], [n,n])
+
+    call check(all(a == b), &
+          msg="all(a == b) failed.",warn=warn)
+
+    call check(sum(diag(a,-1)) - (n-1) < epsilon(1.0_sp), &
+          msg="sum(diag(a,-1)) - (n-1) < epsilon(1.0_sp) failed.",warn=warn)
+
+    call check(all(a == transpose(diag([(1._sp,i=1,n-1)],1))), &
+          msg="all(a == transpose(diag([(1._sp,i=1,n-1)],1))) failed",warn=warn)
+
+    call random_number(a)
+    do i = 1, n
+      call check(size(diag(a,i)) == n-i, &
+            msg="size(diag(a,i)) == n-i failed.",warn=warn)
+    end do
+    call check(size(diag(a,n+1)) == 0, &
+        msg="size(diag(a,n+1)) == 0 failed.",warn=warn)
+  end subroutine
+
+  subroutine test_diag_rdp
+    integer, parameter :: n = 3
+    real(dp) :: v(n), a(n,n), b(n,n)
+    integer :: i,j
+    write(*,*) "test_diag_rdp"
+    v = [(i,i=1,n)]
+    a = diag(v)
+    b = reshape([((merge(i,0,i==j), i=1,n), j=1,n)], [n,n])
+    call check(all(a == b), &
+            msg="all(a == b) failed.",warn=warn)
+
+    call check(all(diag(3*a) == 3*v), &
+            msg="all(diag(3*a) == 3*v) failed.",warn=warn)
+  end subroutine
+
+  subroutine test_diag_rqp
+    integer, parameter :: n = 3
+    real(qp) :: v(n), a(n,n), b(n,n)
+    integer :: i,j
+    write(*,*) "test_diag_rqp"
+    v = [(i,i=1,n)]
+    a = diag(v)
+    b = reshape([((merge(i,0,i==j), i=1,n), j=1,n)], [n,n])
+    call check(all(a == b), &
+            msg="all(a == b) failed.", warn=warn)
+
+    call check(all(diag(3*a) == 3*v), &
+            msg="all(diag(3*a) == 3*v) failed.", warn=warn)
+  end subroutine
+
+  subroutine test_diag_csp
+    integer, parameter :: n = 3
+    complex(sp) :: v(n), a(n,n), b(n,n)
+    complex(sp), parameter :: i_ = complex(0,1)
+    integer :: i,j
+    write(*,*) "test_diag_csp"
+    a = diag([(i,i=1,n)]) + diag([(i_,i=1,n)])
+    b = reshape([((merge(i + 1*i_,0*i_,i==j), i=1,n), j=1,n)], [n,n])
+    call check(all(a == b), &
+            msg="all(a == b) failed.",warn=warn)
+
+    call check(all(abs(real(diag(a)) - [(i,i=1,n)]) < epsilon(1.0_sp)), &
+          msg="all(abs(real(diag(a)) - [(i,i=1,n)]) < epsilon(1.0_sp))", warn=warn)
+    call check(all(abs(aimag(diag(a)) - [(1,i=1,n)]) < epsilon(1.0_sp)), &
+          msg="all(abs(aimag(diag(a)) - [(1,i=1,n)]) < epsilon(1.0_sp))", warn=warn)
+  end subroutine
+
+  subroutine test_diag_cdp
+    integer, parameter :: n = 3
+    complex(dp) :: v(n), a(n,n), b(n,n)
+    complex(dp), parameter :: i_ = complex(0,1)
+    integer :: i,j
+    write(*,*) "test_diag_cdp"
+    a = diag([i_],-2) + diag([i_],2)
+    call check(a(3,1) == i_ .and. a(1,3) == i_, &
+          msg="a(3,1) == i_ .and. a(1,3) == i_ failed.",warn=warn)
+  end subroutine
+
+  subroutine test_diag_cqp
+    integer, parameter :: n = 3
+    complex(qp) :: v(n), a(n,n), b(n,n)
+    complex(qp), parameter :: i_ = complex(0,1)
+    integer :: i,j
+    write(*,*) "test_diag_cqp"
+    a = diag([i_,i_],-1) + diag([i_,i_],1)
+    call check(all(diag(a,-1) == i_) .and. all(diag(a,1) == i_), &
+          msg="all(diag(a,-1) == i_) .and. all(diag(a,1) == i_) failed.",warn=warn)
+  end subroutine
+
+  subroutine test_diag_int8
+    integer, parameter :: n = 3
+    integer(int8), allocatable :: a(:,:)
+    integer :: i
+    logical, allocatable :: mask(:,:)
+    write(*,*) "test_diag_int8"
+    a = reshape([(i,i=1,n**2)],[n,n])
+    mask = merge(.true.,.false.,eye(n) == 1)
+    call check(all(diag(a) == pack(a,mask)), &
+          msg="all(diag(a) == pack(a,mask)) failed.", warn=warn)
+    call check(all(diag(diag(a)) == merge(a,0_int8,mask)), &
+          msg="all(diag(diag(a)) == merge(a,0_int8,mask)) failed.", warn=warn)
+  end subroutine
+  subroutine test_diag_int16
+    integer, parameter :: n = 4
+    integer(int16), allocatable :: a(:,:)
+    integer :: i
+    logical, allocatable :: mask(:,:)
+    write(*,*) "test_diag_int16"
+    a = reshape([(i,i=1,n**2)],[n,n])
+    mask = merge(.true.,.false.,eye(n) == 1)
+    call check(all(diag(a) == pack(a,mask)), &
+          msg="all(diag(a) == pack(a,mask))", warn=warn)
+    call check(all(diag(diag(a)) == merge(a,0_int16,mask)), &
+          msg="all(diag(diag(a)) == merge(a,0_int16,mask)) failed.", warn=warn)
+    a = unpack(int([1,2,3,4],int16),eye(n)==1,a)
+  end subroutine
+  subroutine test_diag_int32
+    integer, parameter :: n = 3
+    integer(int32) :: a(n,n)
+    logical :: mask(n,n)
+    integer :: i, j
+    write(*,*) "test_diag_int32"
+    mask = reshape([((merge(.true.,.false.,i==j+1), i=1,n), j=1,n)], [n,n])
+    a = 0
+    a = unpack([1_int32,1_int32],mask,a)
+    call check(all(diag([1,1],-1) == a), &
+          msg="all(diag([1,1],-1) == a) failed.", warn=warn)
+    call check(all(diag([1,1],1) == transpose(a)), &
+          msg="all(diag([1,1],1) == transpose(a)) failed.", warn=warn)
+  end subroutine
+  subroutine test_diag_int64
+    integer, parameter :: n = 4
+    integer(int64) :: a(n,n), c(2*n-1)
+    logical :: mask(n,n)
+    integer :: i, j
+    
+    write(*,*) "test_diag_int64"
+
+    mask = reshape([((merge(.true.,.false.,i+1==j), i=1,n), j=1,n)], [n,n])
+    a = 0
+    a = unpack([1_int64,1_int64,1_int64],mask,a)
+
+    call check(all(diag([1,1,1],1) == a), &
+          msg="all(diag([1,1,1],1) == a) failed.", warn=warn)
+    call check(all(diag([1,1,1],-1) == transpose(a)), &
+          msg="all(diag([1,1,1],-1) == transpose(a)) failed.", warn=warn)
+
+
+    ! Fill array c with Catalan numbers
+    do i = 0, 2*n-1
+      c(i) = catalan_number(i)
+    end do
+    ! Symmetric Hankel matrix filled with Catalan numbers (det(H) = 1)
+    do i = 1, n
+      do j = 1, n
+        a(i,j) = c(i-1 + (j-1))
+      end do
+    end do
+    call check(all(diag(a,-2) == diag(a,2)), &
+          msg="all(diag(a,-2) == diag(a,2))", warn=warn)
+  end subroutine
+
+
+
+
   subroutine test_trace_rsp
     integer, parameter :: n = 5
     real(sp) :: a(n,n)
+    integer :: i
     write(*,*) "test_trace_rsp"
     a = reshape([(i,i=1,n**2)],[n,n])
     call check(abs(trace(a) - sum(diag(a))) < epsilon(1.0_sp), &
@@ -81,22 +268,24 @@ subroutine test_trace_rsp
   subroutine test_trace_rsp_nonsquare
     integer, parameter :: n = 4
     real(sp) :: a(n,n+1), ans
+    integer :: i
     write(*,*) "test_trace_rsp_nonsquare"
-    a = reshape([(i,i=1,n*(n+1))],[n,n+1])
 
     ! 1   5   9  13  17
     ! 2   6  10  14  18
     ! 3   7  11  15  19
     ! 4   8  12  16  20
-
+    a = reshape([(i,i=1,n*(n+1))],[n,n+1])
     ans = sum([1._sp,6._sp,11._sp,16._sp])
+
     call check(abs(trace(a) - ans) < epsilon(1.0_sp), &
           msg="abs(trace(a) - ans) < epsilon(1.0_sp) failed.",warn=warn)
   end subroutine
 
   subroutine test_trace_rdp
     integer, parameter :: n = 4
     real(dp) :: a(n,n)
+    integer :: i
     write(*,*) "test_trace_rdp"
     a = reshape([(i,i=1,n**2)],[n,n])
     call check(abs(trace(a) - sum(diag(a))) < epsilon(1.0_dp), &
@@ -106,22 +295,24 @@ subroutine test_trace_rdp
   subroutine test_trace_rdp_nonsquare
     integer, parameter :: n = 4
     real(dp) :: a(n,n-1), ans
+    integer :: i
     write(*,*) "test_trace_rdp_nonsquare"
-    a = reshape([(i**2,i=1,n*(n-1))],[n,n-1])
 
     !  1   25   81
     !  4   36  100  
     !  9   49  121  
     ! 16   64  144  
-
+    a = reshape([(i**2,i=1,n*(n-1))],[n,n-1])
     ans = sum([1._dp,36._dp,121._dp])
+
     call check(abs(trace(a) - ans) < epsilon(1.0_dp), &
           msg="abs(trace(a) - ans) < epsilon(1.0_sp) failed.",warn=warn)
   end subroutine
 
   subroutine test_trace_rqp
     integer, parameter :: n = 3
     real(qp) :: a(n,n)
+    integer :: i
     write(*,*) "test_trace_rqp"
     a = reshape([(i,i=1,n**2)],[n,n])
     call check(abs(trace(a) - sum(diag(a))) < epsilon(1.0_qp), &
@@ -205,9 +396,9 @@ subroutine test_trace_int32
   end subroutine
 
   subroutine test_trace_int64
-    integer(int64), parameter :: n = 5
-    integer(int64), parameter :: nd = 2*n-1 ! number of diagonals
-    integer(int64) :: i, j
+    integer, parameter :: n = 5
+    integer, parameter :: nd = 2*n-1 ! number of diagonals
+    integer :: i, j
     integer(int64) :: c(0:nd), H(n,n)
     write(*,*) "test_trace_int64"
 
@@ -229,15 +420,15 @@ subroutine test_trace_int64
   end subroutine
 
   pure recursive function catalan_number(n) result(value)
-    integer(int64), intent(in) :: n
-    integer(int64) :: value
+    integer, intent(in) :: n
+    integer :: value
     integer :: i
     if (n <= 1) then
       value = 1
     else
       value = 0
       do i = 0, n-1
-        value = value + catalan_number(int(i,int64))*catalan_number(n-i-1)
+        value = value + catalan_number(i)*catalan_number(n-i-1)
       end do
     end if
   end function