|
| 1 | +#:include "common.fypp" |
| 2 | +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES |
| 3 | +module test_linalg_mnorm |
| 4 | + use testdrive, only: error_type, check, new_unittest, unittest_type |
| 5 | + use stdlib_linalg_constants |
| 6 | + use stdlib_linalg, only: mnorm, linalg_state_type |
| 7 | + use stdlib_linalg_state, only: linalg_state_type |
| 8 | + |
| 9 | + implicit none (type,external) |
| 10 | + |
| 11 | + contains |
| 12 | + |
| 13 | + !> Matrix norm tests |
| 14 | + subroutine test_matrix_norms(tests) |
| 15 | + !> Collection of tests |
| 16 | + type(unittest_type), allocatable, intent(out) :: tests(:) |
| 17 | + |
| 18 | + allocate(tests(0)) |
| 19 | + |
| 20 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 21 | + tests = [tests,new_unittest("test_matrix_norms_${ri}$",test_matrix_norms_${ri}$)] |
| 22 | + #:endfor |
| 23 | + |
| 24 | + end subroutine test_matrix_norms |
| 25 | + |
| 26 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 27 | + !> Test 1-norm, 2-norm (Euclidean), and infinity norm for ${rt}$ matrices |
| 28 | + subroutine test_matrix_norms_${ri}$(error) |
| 29 | + type(error_type), allocatable, intent(out) :: error |
| 30 | + |
| 31 | + integer(ilp) :: i,j,m,n |
| 32 | + integer(ilp), parameter :: mtx_dim = 5 |
| 33 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 34 | + ${rt}$, allocatable :: A(:,:) |
| 35 | + type(linalg_state_type) :: err |
| 36 | + character(64) :: msg |
| 37 | + |
| 38 | + allocate(A(mtx_dim,mtx_dim)) |
| 39 | + |
| 40 | + ! Initialize matrix with small values to avoid overflow |
| 41 | + A = reshape([(0.01_${rk}$*(i-mtx_dim/2_ilp), i=1_ilp,mtx_dim*mtx_dim)], [mtx_dim,mtx_dim]) |
| 42 | + |
| 43 | + ! 1-norm (Maximum absolute column sum) |
| 44 | + call check(error, abs(mnorm(A, '1', err) - maxval(sum(abs(A), dim=1),1)) < tol*mnorm(A, '1', err), & |
| 45 | + 'Matrix 1-norm does not match expected value') |
| 46 | + if (allocated(error)) return |
| 47 | + |
| 48 | + ! 2-norm (Frobenius norm) |
| 49 | + call check(error, abs(mnorm(A, '2', err) - sqrt(sum(A**2)) < tol*mnorm(A, '2', err), & |
| 50 | + 'Matrix Frobenius norm does not match expected value') |
| 51 | + if (allocated(error)) return |
| 52 | + |
| 53 | + ! Inf-norm (Maximum absolute row sum) |
| 54 | + call check(error, abs(mnorm(A, 'Inf', err) - maxval(sum(abs(A), dim=2),1)) < tol*mnorm(A, 'Inf', err), & |
| 55 | + 'Matrix Infinity norm does not match expected value') |
| 56 | + if (allocated(error)) return |
| 57 | + |
| 58 | + end subroutine test_matrix_norms_${ri}$ |
| 59 | + |
| 60 | + #:for rank in range(3, MAXRANK) |
| 61 | + !> Test N-D norms |
| 62 | + subroutine test_mnorm_${ri}$_${rank}$d(error) |
| 63 | + type(error_type), allocatable, intent(out) :: error |
| 64 | + |
| 65 | + integer(ilp) :: j,dim(2) |
| 66 | + integer(ilp), parameter :: ndim = ${rank}$ |
| 67 | + integer(ilp), parameter :: n = 2_ilp**ndim |
| 68 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 69 | + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ |
| 70 | + |
| 71 | + character(64) :: msg |
| 72 | + |
| 73 | + allocate(a(n), b${fixedranksuffix(rank,2)}$) |
| 74 | + |
| 75 | + ! Init as a range,but with small elements such that all power norms will |
| 76 | + ! never overflow, even in single precision |
| 77 | + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] |
| 78 | + b = reshape(a, shape(b)) |
| 79 | + |
| 80 | + |
| 81 | + |
| 82 | + ! Test norm as collapsed around dimension |
| 83 | + do dim = 1, ndim |
| 84 | + write(msg,"('Not all dim=',i0,' Euclidean norms match ${rt}$ `norm2` intrinsic')") dim |
| 85 | + call check(error,all(abs(mnorm(b,2,dim)-norm2(b,dim))<tol*max(1.0_${rk}$,norm(b,2,dim))),& |
| 86 | + trim(msg)) |
| 87 | + if (allocated(error)) return |
| 88 | + end do |
| 89 | + |
| 90 | + end subroutine test_mnorm_${ri}$_${rank}$d |
| 91 | + #:endif |
| 92 | + |
| 93 | + |
| 94 | + #:endfor |
| 95 | + |
| 96 | +end module test_linalg_mnorm |
| 97 | + |
| 98 | +program test_mnorm |
| 99 | + use, intrinsic :: iso_fortran_env, only : error_unit |
| 100 | + use testdrive, only : run_testsuite, new_testsuite, testsuite_type |
| 101 | + use test_linalg_mnorm, only : test_matrix_norms |
| 102 | + implicit none |
| 103 | + integer :: stat, is |
| 104 | + type(testsuite_type), allocatable :: testsuites(:) |
| 105 | + character(len=*), parameter :: fmt = '("#", *(1x, a))' |
| 106 | + |
| 107 | + stat = 0 |
| 108 | + |
| 109 | + testsuites = [ & |
| 110 | + new_testsuite("matrix_norms", test_matrix_norms) & |
| 111 | + ] |
| 112 | + |
| 113 | + do is = 1, size(testsuites) |
| 114 | + write(error_unit, fmt) "Testing:", testsuites(is)%name |
| 115 | + call run_testsuite(testsuites(is)%collect, error_unit, stat) |
| 116 | + end do |
| 117 | + |
| 118 | + if (stat > 0) then |
| 119 | + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" |
| 120 | + error stop |
| 121 | + end if |
| 122 | +end program test_mnorm |
0 commit comments