|
| 1 | +#:include "common.fypp" |
| 2 | +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES |
| 3 | +! Test Cholesky factorization |
| 4 | +module test_linalg_cholesky |
| 5 | + use testdrive, only: error_type, check, new_unittest, unittest_type |
| 6 | + use stdlib_linalg_constants |
| 7 | + use stdlib_linalg, only: cholesky |
| 8 | + use stdlib_linalg_state, only: linalg_state_type |
| 9 | + |
| 10 | + implicit none (type,external) |
| 11 | + private |
| 12 | + |
| 13 | + public :: test_cholesky_factorization |
| 14 | + |
| 15 | + contains |
| 16 | + |
| 17 | + !> Cholesky factorization tests |
| 18 | + subroutine test_cholesky_factorization(tests) |
| 19 | + !> Collection of tests |
| 20 | + type(unittest_type), allocatable, intent(out) :: tests(:) |
| 21 | + |
| 22 | + allocate(tests(0)) |
| 23 | + |
| 24 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 25 | + #:if rk!="xdp" |
| 26 | + tests = [tests,new_unittest("least_cholesky_${ri}$",test_cholesky_${ri}$)] |
| 27 | + #:endif |
| 28 | + #:endfor |
| 29 | + |
| 30 | + end subroutine test_cholesky_factorization |
| 31 | + |
| 32 | + !> Cholesky factorization of a random matrix |
| 33 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 34 | + #:if rk!="xdp" |
| 35 | + subroutine test_cholesky_${ri}$(error) |
| 36 | + type(error_type), allocatable, intent(out) :: error |
| 37 | + |
| 38 | + integer(ilp), parameter :: n = 3_ilp |
| 39 | + real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$)) |
| 40 | + ${rt}$ :: a(n,n), l(n,n) |
| 41 | + type(linalg_state_type) :: state |
| 42 | + |
| 43 | + ! Set real matrix |
| 44 | + a(1,:) = [6, 15, 55] |
| 45 | + a(2,:) = [15, 55, 225] |
| 46 | + a(3,:) = [55, 225, 979] |
| 47 | + |
| 48 | + ! Set result (lower factor) |
| 49 | + l(1,:) = [ 2.4495_${rk}$, 0.0000_${rk}$, 0.0000_${rk}$] |
| 50 | + l(2,:) = [ 6.1237_${rk}$, 4.1833_${rk}$, 0.0000_${rk}$] |
| 51 | + l(3,:) = [22.4537_${rk}$, 20.9165_${rk}$, 6.1101_${rk}$] |
| 52 | + |
| 53 | + ! 1) Cholesky factorization with full matrices |
| 54 | + call cholesky(a, l, other_zeroed=.true., err=state) |
| 55 | + |
| 56 | + call check(error,state%ok(),state%print()) |
| 57 | + if (allocated(error)) return |
| 58 | + |
| 59 | + call check(error, all(abs(a-matmul(l,transpose(l)))<tol), 'data converged') |
| 60 | + if (allocated(error)) return |
| 61 | + |
| 62 | + |
| 63 | + end subroutine test_cholesky_${ri}$ |
| 64 | + |
| 65 | + #:endif |
| 66 | + #:endfor |
| 67 | + |
| 68 | + |
| 69 | +end module test_linalg_cholesky |
| 70 | + |
| 71 | +program test_cholesky |
| 72 | + use, intrinsic :: iso_fortran_env, only : error_unit |
| 73 | + use testdrive, only : run_testsuite, new_testsuite, testsuite_type |
| 74 | + use test_linalg_cholesky, only : test_cholesky_factorization |
| 75 | + implicit none |
| 76 | + integer :: stat, is |
| 77 | + type(testsuite_type), allocatable :: testsuites(:) |
| 78 | + character(len=*), parameter :: fmt = '("#", *(1x, a))' |
| 79 | + |
| 80 | + stat = 0 |
| 81 | + |
| 82 | + testsuites = [ & |
| 83 | + new_testsuite("linalg_cholesky", test_cholesky_factorization) & |
| 84 | + ] |
| 85 | + |
| 86 | + do is = 1, size(testsuites) |
| 87 | + write(error_unit, fmt) "Testing:", testsuites(is)%name |
| 88 | + call run_testsuite(testsuites(is)%collect, error_unit, stat) |
| 89 | + end do |
| 90 | + |
| 91 | + if (stat > 0) then |
| 92 | + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" |
| 93 | + error stop |
| 94 | + end if |
| 95 | +end program test_cholesky |
0 commit comments