|
| 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 | + #:for rank in range(3, MAXRANK) |
| 23 | + tests = [tests,new_unittest("test_mnorm_${ri}$_${rank}$d",test_mnorm_${ri}$_${rank}$d)] |
| 24 | + #:endfor |
| 25 | + #:endfor |
| 26 | + |
| 27 | + end subroutine test_matrix_norms |
| 28 | + |
| 29 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 30 | + !> Test 1-norm, 2-norm (Euclidean), and infinity norm for ${rt}$ matrices |
| 31 | + subroutine test_matrix_norms_${ri}$(error) |
| 32 | + type(error_type), allocatable, intent(out) :: error |
| 33 | + |
| 34 | + integer(ilp) :: i,j,m,n |
| 35 | + integer(ilp), parameter :: mtx_dim = 5 |
| 36 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 37 | + ${rt}$, allocatable :: A(:,:) |
| 38 | + type(linalg_state_type) :: err |
| 39 | + character(64) :: msg |
| 40 | + |
| 41 | + allocate(A(mtx_dim,mtx_dim)) |
| 42 | + |
| 43 | + ! Initialize matrix with small values to avoid overflow |
| 44 | + A = reshape([(0.01_${rk}$*(i-mtx_dim/2_ilp), i=1_ilp,mtx_dim*mtx_dim)], [mtx_dim,mtx_dim]) |
| 45 | + |
| 46 | + ! 1-norm (Maximum absolute column sum) |
| 47 | + call check(error, abs(mnorm(A, '1', err) - maxval(sum(abs(A), dim=1),1)) < tol*mnorm(A, '1', err), & |
| 48 | + 'Matrix 1-norm does not match expected value') |
| 49 | + if (allocated(error)) return |
| 50 | + |
| 51 | + ! 2-norm (Frobenius norm) |
| 52 | + call check(error, abs(mnorm(A, '2', err) - sqrt(sum(A**2))) < tol*mnorm(A, '2', err), & |
| 53 | + 'Matrix Frobenius norm does not match expected value') |
| 54 | + if (allocated(error)) return |
| 55 | + |
| 56 | + ! Inf-norm (Maximum absolute row sum) |
| 57 | + call check(error, abs(mnorm(A, 'Inf', err) - maxval(sum(abs(A), dim=2),1)) < tol*mnorm(A, 'Inf', err), & |
| 58 | + 'Matrix Infinity norm does not match expected value') |
| 59 | + if (allocated(error)) return |
| 60 | + |
| 61 | + end subroutine test_matrix_norms_${ri}$ |
| 62 | + |
| 63 | + #:for rank in range(3, MAXRANK) |
| 64 | + !> Test N-D norms |
| 65 | + subroutine test_mnorm_${ri}$_${rank}$d(error) |
| 66 | + type(error_type), allocatable, intent(out) :: error |
| 67 | + |
| 68 | + integer(ilp) :: j,dim1,dim2,dim(2),order |
| 69 | + integer(ilp), parameter :: orders(*) = [1_ilp,2_ilp,huge(0_ilp)] |
| 70 | + integer(ilp), parameter :: ndim = ${rank}$ |
| 71 | + integer(ilp), parameter :: n = 2_ilp**ndim |
| 72 | + integer(ilp), parameter :: dims(*) = [(dim1, dim1=1,ndim)] |
| 73 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 74 | + real(${rk}$), allocatable :: bnrm${ranksuffix(rank-2)}$ |
| 75 | + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ |
| 76 | + |
| 77 | + character(64) :: msg |
| 78 | + |
| 79 | + allocate(a(n), b${fixedranksuffix(rank,2)}$) |
| 80 | + |
| 81 | + ! Init as a range,but with small elements such that all power norms will |
| 82 | + ! never overflow, even in single precision |
| 83 | + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] |
| 84 | + b = reshape(a, shape(b)) |
| 85 | + |
| 86 | + ! Test norm as collapsed around dimensions |
| 87 | + do j = 1, size(orders) |
| 88 | + order = orders(j) |
| 89 | + do dim1 = 1, ndim |
| 90 | + do dim2 = dim1+1, ndim |
| 91 | + |
| 92 | + dim = [dim1,dim2] |
| 93 | + |
| 94 | + ! Get norms collapsed on these dims |
| 95 | + bnrm = mnorm(b,order,dim) |
| 96 | + |
| 97 | + ! Assert size |
| 98 | + write(msg,"('dim=[',i0,',',i0,'] order=',i0,' ${rk}$ norm returned wrong shape')") dim, order |
| 99 | + call check(error,all(shape(bnrm)==pack(shape(b),dims/=dim1 .and. dims/=dim2) ), trim(msg)) |
| 100 | + if (allocated(error)) return |
| 101 | + |
| 102 | + end do |
| 103 | + end do |
| 104 | + end do |
| 105 | + |
| 106 | + end subroutine test_mnorm_${ri}$_${rank}$d |
| 107 | + #:endfor |
| 108 | + |
| 109 | + |
| 110 | + #:endfor |
| 111 | + |
| 112 | +end module test_linalg_mnorm |
| 113 | + |
| 114 | +program test_mnorm |
| 115 | + use, intrinsic :: iso_fortran_env, only : error_unit |
| 116 | + use testdrive, only : run_testsuite, new_testsuite, testsuite_type |
| 117 | + use test_linalg_mnorm, only : test_matrix_norms |
| 118 | + implicit none |
| 119 | + integer :: stat, is |
| 120 | + type(testsuite_type), allocatable :: testsuites(:) |
| 121 | + character(len=*), parameter :: fmt = '("#", *(1x, a))' |
| 122 | + |
| 123 | + stat = 0 |
| 124 | + |
| 125 | + testsuites = [ & |
| 126 | + new_testsuite("matrix_norms", test_matrix_norms) & |
| 127 | + ] |
| 128 | + |
| 129 | + do is = 1, size(testsuites) |
| 130 | + write(error_unit, fmt) "Testing:", testsuites(is)%name |
| 131 | + call run_testsuite(testsuites(is)%collect, error_unit, stat) |
| 132 | + end do |
| 133 | + |
| 134 | + if (stat > 0) then |
| 135 | + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" |
| 136 | + error stop |
| 137 | + end if |
| 138 | +end program test_mnorm |
0 commit comments