new test: random SPD matrix (real, complex)

Author: perazz

test/linalg/test_linalg_inverse.fypp

@@ -24,7 +24,8 @@ module test_linalg_inverse
         #:for rk,rt,ri in RC_KINDS_TYPES
         #:if rk!="xdp"
         tests = [tests,new_unittest("${ri}$_eye_inverse",test_${ri}$_eye_inverse), &
-                       new_unittest("${ri}$_singular_inverse",test_${ri}$_singular_inverse)]
+                       new_unittest("${ri}$_singular_inverse",test_${ri}$_singular_inverse), &
+                       new_unittest("${ri}$_random_spd_inverse",test_${ri}$_random_spd_inverse)]
         #:endif
         #:endfor
 
@@ -91,7 +92,53 @@ module test_linalg_inverse
         if (allocated(error)) return
 
     end subroutine test_${ri}$_singular_inverse
+    
+    !> Create a random symmetric positive definite matrix
+    function random_spd_matrix_${ri}$(n) result(A)
+        integer(ilp), intent(in) :: n
+        ${rt}$ :: A(n,n)
+        
+        ${rt}$, parameter :: one  = 1.0_${rk}$
+        ${rt}$, parameter :: half = 0.5_${rk}$
+        
+        !> Initialize with randoms
+        call random_number(A)
+        
+        !> Make symmetric
+        A = half*(A+transpose(A))
+        
+        !> Add diagonally dominant part 
+        A = A + n*eye(n)
+        
+    end function random_spd_matrix_${ri}$
 
+    !> Test random symmetric positive-definite matrix
+    subroutine test_${ri}$_random_spd_inverse(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        !> Solution tolerance
+        ${rt}$, parameter :: tol = sqrt(epsilon(0.0_${rk}$))
+
+        !> Local variables
+        integer(ilp), parameter :: n = 5_ilp
+        type(linalg_state_type) :: state
+        ${rt}$ :: A(n,n),Am1(n,n)
+        
+        !> Generate random SPD matrix
+        A = random_spd_matrix_${ri}$(n)
+
+        !> Invert matrix
+        call invert(A,Am1,err=state)
+        
+        !> Check result    
+        call check(error,state%ok(),'random SPD matrix (${rk}$): '//state%print())
+        if (allocated(error)) return     
+
+        call check(error,all(abs(matmul(Am1,A)-eye(n))<tol),'random SPD matrix (${rk}$): accuracy test')
+        if (allocated(error)) return    
+    
+    end subroutine test_${ri}$_random_spd_inverse
+    
     #:endif
     #:endfor
 
@@ -160,6 +207,76 @@ module test_linalg_inverse
 
     end subroutine test_${ci}$_eye_inverse
 
+    !> Create a random symmetric positive definite matrix
+    function random_spd_matrix_${ci}$(n) result(A)
+        integer(ilp), intent(in) :: n
+        ${ct}$ :: A(n,n)
+        
+        ${ct}$, parameter :: half = (0.5_${ck}$,0.0_${ck}$)
+        real(${ck}$) :: reA(n,n),imA(n,n)
+        integer(ilp) :: i
+        
+        !> Initialize with randoms
+        call random_number(reA)
+        call random_number(imA)
+        
+        A = cmplx(reA,imA,kind=${ck}$)
+        
+        !> Make symmetric
+        A = half*(A+transpose(A))
+        
+        !> Add diagonally dominant part 
+        forall(i=1:n) A(i,i) = A(i,i) + n*(1.0_${ck}$,0.0_${ck}$)
+        
+    end function random_spd_matrix_${ci}$
+
+    !> Test random symmetric positive-definite matrix
+    subroutine test_${ci}$_random_spd_inverse(error)
+        type(error_type), allocatable, intent(out) :: error
+
+        !> Local variables
+        integer(ilp) :: failed,i,j
+        integer(ilp), parameter :: n = 5_ilp
+        type(linalg_state_type) :: state
+        ${ct}$ :: A(n,n),Am1(n,n),AA(n,n)
+        
+        !> Generate random SPD matrix
+        A = random_spd_matrix_${ci}$(n)
+
+        !> Invert matrix
+        call invert(A,Am1,err=state)
+        
+        !> Check result    
+        call check(error,state%ok(),'random complex SPD matrix (${ck}$): '//state%print())
+        if (allocated(error)) return     
+
+        failed = 0
+        AA = matmul(A,Am1)
+        do i=1,n
+            do j=1,n
+                if (.not.is_complex_inverse(AA(i,j),i,j)) failed = failed+1
+            end do
+        end do
+
+        call check(error,failed==0,'inverse_${ci}$_eye (subroutine): data converged')
+        if (allocated(error)) return   
+
+        contains
+
+           elemental logical function is_complex_inverse(aij,i,j)
+               ${ct}$, intent(in) :: aij
+               integer(ilp), intent(in) :: i,j
+               real(${ck}$), parameter :: tol = sqrt(epsilon(0.0_${ck}$))
+               if (i/=j) then
+                  is_complex_inverse = abs(aij)<tol
+               else
+                  ! Product should return the real identity
+                  is_complex_inverse = abs(aij - (1.0_${ck}$,0.0_${ck}$))<tol
+               endif
+           end function is_complex_inverse
+    
+    end subroutine test_${ci}$_random_spd_inverse
+
     !> Invert singular matrix
     subroutine test_${ci}$_singular_inverse(error)
         type(error_type), allocatable, intent(out) :: error