real and complex tests

test/linalg/CMakeLists.txt

@@ -2,10 +2,12 @@ set(
  fppFiles
  "test_linalg.fypp"
  "test_blas_lapack.fypp"
+ "test_linalg_solve.fypp"
  "test_linalg_matrix_property_checks.fypp"
 )
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
 
 ADDTEST(linalg)
 ADDTEST(linalg_matrix_property_checks)
+ADDTEST(linalg_solve)
 ADDTEST(blas_lapack)

test/linalg/test_linalg_solve.fypp

@@ -0,0 +1,190 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+! Test linear system solver
+module test_linalg_solve
+ use stdlib_linalg_constants
+ use stdlib_linalg_state
+ use stdlib_linalg, only: solve
+ use testdrive, only: error_type, check, new_unittest, unittest_type 
+
+ implicit none (type,external)
+ private
+
+ public :: test_linear_systems
+
+ contains
+
+ !> Solve real and complex linear systems
+ subroutine test_linear_systems(tests)
+ !> Collection of tests
+ type(unittest_type), allocatable, intent(out) :: tests(:)
+
+ allocate(tests(0))
+
+ #:for rk,rt,ri in REAL_KINDS_TYPES
+ #:if rk!="xdp"
+ tests = [tests,new_unittest("solve_${ri}$",test_${ri}$_solve), &
+ new_unittest("solve_${ri}$_multiple",test_${ri}$_solve_multiple)]
+ #:endif
+ #:endfor
+
+ #:for ck,ct,ci in CMPLX_KINDS_TYPES
+ #:if ck!="xdp"
+ tests = [tests,new_unittest("solve_complex_${ci}$",test_${ci}$_solve), &
+ new_unittest("solve_2x2_complex_${ci}$",test_2x2_${ci}$_solve)]
+ #:endif
+ #:endfor
+
+ end subroutine test_linear_systems 
+
+ #:for rk,rt,ri in REAL_KINDS_TYPES
+ #:if rk!="xdp" 
+ !> Simple linear system
+ subroutine test_${ri}$_solve(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ type(linalg_state_type) :: state
+
+ ${rt}$ :: A(3,3) = transpose(reshape([${rt}$ :: 1, 3, 3, &
+ 1, 3, 4, &
+ 1, 4, 3], [3,3]))
+ ${rt}$ :: b (3) = [${rt}$ :: 1, 4, -1]
+ ${rt}$ :: res(3) = [${rt}$ :: -2, -2, 3]
+ ${rt}$ :: x(3)
+
+ x = solve(a,b,err=state)
+
+ call check(error,state%ok(),state%print())
+ if (allocated(error)) return
+
+ call check(error, all(abs(x-res)<abs(res*epsilon(0.0_${rk}$))), 'results match expected')
+ if (allocated(error)) return
+
+ end subroutine test_${ri}$_solve
+
+ !> Simple linear system with multiple right hand sides
+ subroutine test_${ri}$_solve_multiple(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ type(linalg_state_type) :: state
+
+ ${rt}$ :: A(3,3) = transpose(reshape([${rt}$ :: 1,-1, 2, &
+ 0, 1, 1, &
+ 1,-1, 3], [3,3]))
+ ${rt}$ :: b(3,3) = transpose(reshape([${rt}$ :: 0, 1, 2, &
+ 1,-2,-1, &
+ 2, 3,-1], [3,3]))
+ ${rt}$ :: res(3,3) = transpose(reshape([${rt}$ ::-5,-7,10, &
+ -1,-4, 2, &
+ 2, 2,-3], [3,3]))
+ ${rt}$ :: x(3,3)
+
+ x = solve(a,b,err=state)
+
+ call check(error,state%ok(),state%print())
+ if (allocated(error)) return
+
+ call check(error, all(abs(x-res)<abs(res*epsilon(0.0_${rk}$))), 'results match expected')
+ if (allocated(error)) return 
+
+ end subroutine test_${ri}$_solve_multiple
+ #:endif
+ #:endfor
+
+ #:for rk,rt,ri in CMPLX_KINDS_TYPES
+ #:if rk!="xdp"
+ !> Complex linear system
+ !> Militaru, Popa, "On the numerical solving of complex linear systems",
+ !> Int J Pure Appl Math 76(1), 113-122, 2012.
+ subroutine test_${ri}$_solve(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ type(linalg_state_type) :: state
+
+ ${rt}$ :: A(5,5), b(5), res(5), x(5)
+ integer(ilp) :: i
+
+ ! Fill in linear system
+ A = (0.0_${rk}$,0.0_${rk}$)
+
+ A(1:2,1) = [(19.73_${rk}$,0.0_${rk}$),(0.0_${rk}$,-0.51_${rk}$)]
+ A(1:3,2) = [(12.11_${rk}$,-1.0_${rk}$),(32.3_${rk}$,7.0_${rk}$),(0.0_${rk}$,-0.51_${rk}$)]
+ A(1:4,3) = [(0.0_${rk}$,5.0_${rk}$),(23.07_${rk}$,0.0_${rk}$),(70.0_${rk}$,7.3_${rk}$),(1.0_${rk}$,1.1_${rk}$)]
+ A(2:5,4) = [(0.0_${rk}$,1.0_${rk}$),(3.95_${rk}$,0.0_${rk}$),(50.17_${rk}$,0.0_${rk}$),(0.0_${rk}$,-9.351_${rk}$)]
+ A(3:5,5) = [(19.0_${rk}$,31.83_${rk}$),(45.51_${rk}$,0.0_${rk}$),(55.0_${rk}$,0.0_${rk}$)]
+
+ b = [(77.38_${rk}$,8.82_${rk}$),(157.48_${rk}$,19.8_${rk}$),(1175.62_${rk}$,20.69_${rk}$),(912.12_${rk}$,-801.75_${rk}$),(550.0_${rk}$,-1060.4_${rk}$)]
+
+ ! Exact result
+ res = [(3.3_${rk}$,-1.0_${rk}$),(1.0_${rk}$,0.17_${rk}$),(5.5_${rk}$,0.0_${rk}$),(9.0_${rk}$,0.0_${rk}$),(10.0_${rk}$,-17.75_${rk}$)]
+
+ x = solve(a,b,err=state)
+
+ call check(error,state%ok(),state%print())
+ if (allocated(error)) return
+
+ call check(error, all(abs(x-res)<abs(res)*1.0e-3_${rk}$), 'results match expected')
+ if (allocated(error)) return 
+
+ end subroutine test_${ri}$_solve
+
+ !> 2x2 Complex linear system
+ !> https://math.stackexchange.com/questions/1996540/solving-linear-equation-systems-with-complex-coefficients-and-variables
+ subroutine test_2x2_${ri}$_solve(error)
+ type(error_type), allocatable, intent(out) :: error
+
+ type(linalg_state_type) :: state
+
+ ${rt}$ :: A(2,2), b(2), res(2), x(2)
+ integer(ilp) :: i
+
+ ! Fill in linear system
+ A(1,:) = [(+1.0_${rk}$,+1.0_${rk}$),(-1.0_${rk}$,0.0_${rk}$)]
+ A(2,:) = [(+1.0_${rk}$,-1.0_${rk}$),(+1.0_${rk}$,1.0_${rk}$)]
+
+ b = [(0.0_${rk}$,1.0_${rk}$),(1.0_${rk}$,0.0_${rk}$)]
+
+ ! Exact result
+ res = [(0.5_${rk}$,0.5_${rk}$),(0.0_${rk}$,0.0_${rk}$)]
+
+ x = solve(a,b,err=state)
+
+ call check(error,state%ok(),state%print())
+ if (allocated(error)) return
+
+ call check(error, all(abs(x-res)<max(tiny(0.0_${rk}$),abs(res)*epsilon(0.0_${rk}$))), 'results match expected')
+ if (allocated(error)) return 
+
+
+ end subroutine test_2x2_${ri}$_solve
+ #:endif
+ #:endfor
+
+end module test_linalg_solve
+
+program test_solve
+ use, intrinsic :: iso_fortran_env, only : error_unit
+ use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+ use test_linalg_solve, only : test_linear_systems
+ implicit none
+ integer :: stat, is
+ type(testsuite_type), allocatable :: testsuites(:)
+ character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+ stat = 0
+
+ testsuites = [ &
+ new_testsuite("linalg_solve", test_linear_systems) &
+ ]
+
+ do is = 1, size(testsuites)
+ write(error_unit, fmt) "Testing:", testsuites(is)%name
+ call run_testsuite(testsuites(is)%collect, error_unit, stat)
+ end do
+
+ if (stat > 0) then
+ write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+ error stop
+ end if
+end program test_solve
+