|
| 1 | +#:include "common.fypp" |
| 2 | +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES |
| 3 | + |
| 4 | +#! Generate an array rank suffix with the same fixed size for all dimensions. |
| 5 | +#! |
| 6 | +#! Args: |
| 7 | +#! rank (int): Rank of the variable |
| 8 | +#! size (int): Size along each dimension |
| 9 | +#! |
| 10 | +#! Returns: |
| 11 | +#! Array rank suffix string (e.g. (4,4,4) if rank = 3 and size = 4) |
| 12 | +#! |
| 13 | +#:def fixedranksuffix(rank,size) |
| 14 | +#{if rank > 0}#(${str(size) + (","+str(size)) * (rank - 1)}$)#{endif}# |
| 15 | +#:enddef |
| 16 | +! Test vector norms |
| 17 | +module test_linalg_norm |
| 18 | + use testdrive, only: error_type, check, new_unittest, unittest_type |
| 19 | + use stdlib_linalg_constants |
| 20 | + use stdlib_linalg, only: norm, linalg_state_type |
| 21 | + use stdlib_linalg_state, only: linalg_state_type |
| 22 | + |
| 23 | + implicit none (type,external) |
| 24 | + |
| 25 | + contains |
| 26 | + |
| 27 | + !> Vector norm tests |
| 28 | + subroutine test_vector_norms(tests) |
| 29 | + !> Collection of tests |
| 30 | + type(unittest_type), allocatable, intent(out) :: tests(:) |
| 31 | + |
| 32 | + allocate(tests(0)) |
| 33 | + |
| 34 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 35 | + #:for rank in range(1, MAXRANK) |
| 36 | + tests = [tests,new_unittest("norm_${ri}$_${rank}$d",test_norm_${ri}$_${rank}$d)] |
| 37 | + #:endfor |
| 38 | + #:for rank in range(2, MAXRANK) |
| 39 | + #:if rt.startswith('real') |
| 40 | + tests = [tests,new_unittest("norm2_${ri}$_${rank}$d",test_norm2_${ri}$_${rank}$d)] |
| 41 | + #:endif |
| 42 | + tests = [tests,new_unittest("norm_dimmed_${ri}$_${rank}$d",test_norm_dimmed_${ri}$_${rank}$d)] |
| 43 | + #:endfor |
| 44 | + #:endfor |
| 45 | + |
| 46 | + end subroutine test_vector_norms |
| 47 | + |
| 48 | + #:for rk,rt,ri in RC_KINDS_TYPES |
| 49 | + #:for rank in range(1, MAXRANK) |
| 50 | + |
| 51 | + !> Test several norms with different dimensions |
| 52 | + subroutine test_norm_${ri}$_${rank}$d(error) |
| 53 | + type(error_type), allocatable, intent(out) :: error |
| 54 | + |
| 55 | + integer(ilp) :: j,order |
| 56 | + integer(ilp), parameter :: n = 2_ilp**${rank}$ |
| 57 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 58 | + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ |
| 59 | + character(64) :: msg |
| 60 | + |
| 61 | + allocate(a(n), b${fixedranksuffix(rank,2)}$) |
| 62 | + |
| 63 | + ! Init as a range,but with small elements such that all power norms will |
| 64 | + ! never overflow, even in single precision |
| 65 | + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] |
| 66 | + b = reshape(a, shape(b)) |
| 67 | + |
| 68 | + ! Test some norms |
| 69 | + do order = 1, 10 |
| 70 | + write(msg,"('reshaped order-',i0,' p-norm is the same')") order |
| 71 | + call check(error,abs(norm(a,order)-norm(b,order))<tol*max(1.0_${rk}$,norm(a,order)),trim(msg)) |
| 72 | + if (allocated(error)) return |
| 73 | + end do |
| 74 | + |
| 75 | + ! Infinity norms |
| 76 | + call check(error,abs(norm(a,'inf')-norm(b,'inf'))<tol*max(1.0_${rk}$,norm(a,'inf')),& |
| 77 | + 'reshaped +infinity norm is the same') |
| 78 | + if (allocated(error)) return |
| 79 | + |
| 80 | + ! Infinity norms |
| 81 | + call check(error,abs(norm(a,'-inf')-norm(b,'-inf'))<tol*max(1.0_${rk}$,norm(a,'-inf')),& |
| 82 | + 'reshaped -infinity norm is the same') |
| 83 | + if (allocated(error)) return |
| 84 | + |
| 85 | + end subroutine test_norm_${ri}$_${rank}$d |
| 86 | + #:endfor |
| 87 | + |
| 88 | + !> Test Euclidean norm; compare with Fortran intrinsic norm2 for reals |
| 89 | + #:for rank in range(2, MAXRANK) |
| 90 | + #:if rt.startswith('real') |
| 91 | + subroutine test_norm2_${ri}$_${rank}$d(error) |
| 92 | + type(error_type), allocatable, intent(out) :: error |
| 93 | + |
| 94 | + integer(ilp) :: j,dim |
| 95 | + integer(ilp), parameter :: ndim = ${rank}$ |
| 96 | + integer(ilp), parameter :: n = 2_ilp**ndim |
| 97 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 98 | + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ |
| 99 | + intrinsic :: norm2 |
| 100 | + character(64) :: msg |
| 101 | + |
| 102 | + allocate(a(n), b${fixedranksuffix(rank,2)}$) |
| 103 | + |
| 104 | + ! Init as a range,but with small elements such that all power norms will |
| 105 | + ! never overflow, even in single precision |
| 106 | + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] |
| 107 | + b = reshape(a, shape(b)) |
| 108 | + |
| 109 | + ! Test some norms |
| 110 | + call check(error,abs(norm(a,2) - norm2(a))<tol*norm(a,2),& |
| 111 | + 'Euclidean norm does not match `norm2` intrinsic') |
| 112 | + if (allocated(error)) return |
| 113 | + |
| 114 | + ! Infinity norms |
| 115 | + call check(error,abs(norm(b,2)-norm2(b))<tol*norm(b,2),& |
| 116 | + 'Dimmed Euclidean norm does not match `norm2` intrinsic') |
| 117 | + if (allocated(error)) return |
| 118 | + |
| 119 | + ! Test norm as collapsed around dimension |
| 120 | + do dim = 1, ndim |
| 121 | + write(msg,"('Not all dim=',i0,' Euclidean norms match `norm2` intrinsic')") dim |
| 122 | + call check(error,all(abs(norm(b,2,dim)-norm2(b,dim))<tol*max(1.0_${rk}$,norm(b,2,dim))),& |
| 123 | + trim(msg)) |
| 124 | + if (allocated(error)) return |
| 125 | + end do |
| 126 | + |
| 127 | + end subroutine test_norm2_${ri}$_${rank}$d |
| 128 | + #:endif |
| 129 | + |
| 130 | + ! Test norm along a dimension and compare it against individually evaluated norms |
| 131 | + subroutine test_norm_dimmed_${ri}$_${rank}$d(error) |
| 132 | + type(error_type), allocatable, intent(out) :: error |
| 133 | + |
| 134 | + integer(ilp) :: j,dim,order |
| 135 | + integer(ilp), parameter :: ndim = ${rank}$ |
| 136 | + integer(ilp), parameter :: n = 2_ilp**ndim |
| 137 | + integer(ilp), parameter :: dims(*) = [(dim, dim=1,ndim)] |
| 138 | + real(${rk}$), parameter :: tol = 10*sqrt(epsilon(0.0_${rk}$)) |
| 139 | + integer(ilp) :: coords(ndim) |
| 140 | + real :: x(ndim) |
| 141 | + ${rt}$, allocatable :: a(:), b${ranksuffix(rank)}$ |
| 142 | + real(${rk}$), allocatable :: bnrm${ranksuffix(rank-1)}$ |
| 143 | + character(64) :: msg |
| 144 | + |
| 145 | + allocate(a(n), b${fixedranksuffix(rank,2)}$) |
| 146 | + |
| 147 | + ! Init as a range,but with small elements such that all power norms will |
| 148 | + ! never overflow, even in single precision |
| 149 | + a = [(0.01_${rk}$*(j-n/2_ilp), j=1_ilp,n)] |
| 150 | + b = reshape(a, shape(b)) |
| 151 | + |
| 152 | + do order = 1, 5 |
| 153 | + |
| 154 | + do dim = 1, ndim |
| 155 | + |
| 156 | + bnrm = norm(b, order, dim) |
| 157 | + |
| 158 | + ! Assert size |
| 159 | + write(msg,"('dim=',i0,', order=',i0,' ${rk}$ norm has the wrong shape')") dim, order |
| 160 | + call check(error,all( shape(bnrm)==pack(shape(b),dims/=dim) ), trim(msg)) |
| 161 | + if (allocated(error)) return |
| 162 | + |
| 163 | + end do |
| 164 | + |
| 165 | + end do |
| 166 | + |
| 167 | + end subroutine test_norm_dimmed_${ri}$_${rank}$d |
| 168 | + |
| 169 | + |
| 170 | + #:endfor |
| 171 | + #:endfor |
| 172 | + |
| 173 | + |
| 174 | +end module test_linalg_norm |
| 175 | + |
| 176 | +program test_norm |
| 177 | + use, intrinsic :: iso_fortran_env, only : error_unit |
| 178 | + use testdrive, only : run_testsuite, new_testsuite, testsuite_type |
| 179 | + use test_linalg_norm, only : test_vector_norms |
| 180 | + implicit none |
| 181 | + integer :: stat, is |
| 182 | + type(testsuite_type), allocatable :: testsuites(:) |
| 183 | + character(len=*), parameter :: fmt = '("#", *(1x, a))' |
| 184 | + |
| 185 | + stat = 0 |
| 186 | + |
| 187 | + testsuites = [ & |
| 188 | + new_testsuite("linalg_norm", test_vector_norms) & |
| 189 | + ] |
| 190 | + |
| 191 | + do is = 1, size(testsuites) |
| 192 | + write(error_unit, fmt) "Testing:", testsuites(is)%name |
| 193 | + call run_testsuite(testsuites(is)%collect, error_unit, stat) |
| 194 | + end do |
| 195 | + |
| 196 | + if (stat > 0) then |
| 197 | + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" |
| 198 | + error stop |
| 199 | + end if |
| 200 | +end program test_norm |
| 201 | + |
| 202 | + |
0 commit comments