diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 106fa02ac..df4358220 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -6,6 +6,7 @@ set(fppFiles stdlib_experimental_optval.fypp stdlib_experimental_stats.fypp stdlib_experimental_stats_mean.fypp + stdlib_experimental_stats_var.fypp ) diff --git a/src/Makefile.manual b/src/Makefile.manual index f3772f2bd..7159cfa76 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -5,7 +5,8 @@ SRC = stdlib_experimental_ascii.f90 \ stdlib_experimental_kinds.f90 \ f18estop.f90 \ stdlib_experimental_stats.f90 \ - stdlib_experimental_stats_mean.f90 + stdlib_experimental_stats_mean.f90 \ + stdlib_experimental_stats_var.f90 LIB = libstdlib.a @@ -42,6 +43,11 @@ stdlib_experimental_stats_mean.o: \ stdlib_experimental_optval.o \ stdlib_experimental_kinds.o \ stdlib_experimental_stats.o +stdlib_experimental_stats_var.o: \ + stdlib_experimental_optval.o \ + stdlib_experimental_kinds.o \ + stdlib_experimental_stats.o stdlib_experimental_io.f90: stdlib_experimental_io.fypp stdlib_experimental_stats.f90: stdlib_experimental_stats.fypp stdlib_experimental_stats_mean.f90: stdlib_experimental_stats_mean.fypp +stdlib_experimental_stats_var.f90: stdlib_experimental_stats_var.fypp diff --git a/src/common.fypp b/src/common.fypp index 99fed2250..0dbd3e7d0 100644 --- a/src/common.fypp +++ b/src/common.fypp @@ -97,6 +97,7 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ #:endif #:enddef + #! Generates a routine name from a generic name, rank, type and kind #! #! Args: @@ -113,4 +114,28 @@ ${prefix + joinstr.join([line.strip() for line in txt.split("\n")]) + suffix}$ $:"{0}_{1}_{2}{3}_{2}{3}".format(gname, rank, type[0], kind) if suffix == '' else "{0}_{1}_{2}{3}_{4}".format(gname, rank, type[0], kind, suffix) #:enddef + +#! Generates an array rank suffix for subarrays reducing the dimension +#! +#! Args: +#! rank (int): Rank of the original variable +#! selectors (array): Dimension and name of the variable(s) +#! +#! Returns: +#! Array rank suffix string enclosed in braces +#! +#! E.g., +#! select_subarray(5 , [(4, 'i'), (5, 'j')])}$ +#! -> (:, :, :, i, j) +#! +#:def select_subarray(rank, selectors) + #:assert rank > 0 + #:set seldict = dict(selectors) + #:call join_lines(joinstr=", ", prefix="(", suffix=")") + #:for i in range(1, rank + 1) + $:seldict.get(i, ":") + #:endfor + #:endcall +#:enddef + #:endmute diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp index f16fa0a92..ba5e7ad63 100644 --- a/src/stdlib_experimental_stats.fypp +++ b/src/stdlib_experimental_stats.fypp @@ -7,7 +7,7 @@ module stdlib_experimental_stats implicit none private ! Public API - public :: mean + public :: mean, var interface mean #:for k1, t1 in RC_KINDS_TYPES @@ -104,4 +104,102 @@ module stdlib_experimental_stats end interface mean + + interface var + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_all",rank, t1, k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(${k1}$) :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_all",rank, t1, k1, 'dp') + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(dp) :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var",rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var",rank, t1, k1, 'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask_all",rank, t1, k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${k1}$) :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask_all",rank, t1, k1, 'dp') + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask",rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask",rank, t1, k1, 'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + end function ${RName}$ + #:endfor + #:endfor + + end interface var + + end module stdlib_experimental_stats diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md index 0c11c6822..f8fd1fb09 100644 --- a/src/stdlib_experimental_stats.md +++ b/src/stdlib_experimental_stats.md @@ -3,6 +3,7 @@ ## Implemented * `mean` + * `var` ## `mean` - mean of array elements @@ -47,3 +48,55 @@ program demo_mean reshape(x, [ 2, 3 ] ) > 3.) !returns [ NaN, 4.0, 5.5 ] end program demo_mean ``` + +## `var` - variance of array elements + +### Description + +Returns the variance of all the elements of `array`, or of the elements of `array` along dimension `dim` if provided, and if the corresponding element in `mask` is `true`. + +The variance is defined as the best unbiased estimator and is computed as: + +``` + var(x) = 1/(n-1) sum_i (array(i) - mean(array))^2 +``` + +### Syntax + +`result = var(array [, mask])` + +`result = var(array, dim [, mask])` + +### Arguments + +`array`: Shall be an array of type `integer`, `real`, or `complex`. + +`dim`: Shall be a scalar of type `integer` with a value in the range from 1 to n, where n is the rank of `array`. + +`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `array`. + +### Return value + +If `array` is of type `real` or `complex`, the result is of the same type as `array`. +If `array` is of type `integer`, the result is of type `double precision`. + +If `dim` is absent, a scalar with the variance of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `ar +ray` with dimension `dim` dropped is returned. + +If `mask` is specified, the result is the variance of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`. + +### Example + +```fortran +program demo_var + use stdlib_experimental_stats, only: var + implicit none + real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ] + print *, var(x) !returns 3.5 + print *, var( reshape(x, [ 2, 3 ] )) !returns 3.5 + print *, var( reshape(x, [ 2, 3 ] ), 1) !returns [0.5, 0.5, 0.5] + print *, var( reshape(x, [ 2, 3 ] ), 1,& + reshape(x, [ 2, 3 ] ) > 3.) !returns [0., NaN, 0.5] +end program demo_var +``` + diff --git a/src/stdlib_experimental_stats_var.fypp b/src/stdlib_experimental_stats_var.fypp new file mode 100644 index 000000000..eaa5bdcdd --- /dev/null +++ b/src/stdlib_experimental_stats_var.fypp @@ -0,0 +1,265 @@ +#:include "common.fypp" +#:set RANKS = range(1, MAXRANK + 1) +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +submodule (stdlib_experimental_stats) stdlib_experimental_stats_var + + use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan + use stdlib_experimental_error, only: error_stop + use stdlib_experimental_optval, only: optval + implicit none + +contains + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_all",rank, t1, k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(${k1}$) :: res + + real(${k1}$) :: n + ${t1}$ :: mean + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + n = real(size(x, kind = int64), ${k1}$) + mean = sum(x) / n + + #:if t1[0] == 'r' + res = sum((x - mean)**2) / (n - 1._${k1}$) + #:else + res = sum(abs(x - mean)**2) / (n - 1._${k1}$) + #:endif + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_all",rank, t1, k1, 'dp') + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in), optional :: mask + real(dp) :: res + + real(dp) :: n, mean + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + n = real(size(x, kind = int64), dp) + mean = sum(real(x, dp)) / n + + res = sum((real(x, dp) - mean)**2) / (n - 1._dp) + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var",rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(${k1}$) :: n + ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._${k1}$, ieee_quiet_nan) + return + end if + + res = 0._${k1}$ + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(size(x, dim), ${k1}$) + mean = sum(x, dim) / n + do i = 1, size(x, dim) + #:if t1[0] == 'r' + res = res + (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2 + #:else + res = res + abs(x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2 + #:endif + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._${k1}$) + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var",rank, t1, k1, 'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in), optional :: mask + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(dp) :: n + real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + + if (.not.optval(mask, .true.)) then + res = ieee_value(1._dp, ieee_quiet_nan) + return + end if + + res = 0._dp + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(size(x, dim), dp) + mean = sum(real(x, dp), dim) / n + do i = 1, size(x, dim) + res = res + (real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2 + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._dp) + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask_all",rank, t1, k1) + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${k1}$) :: res + + real(${k1}$) :: n + ${t1}$ :: mean + + n = real(count(mask, kind = int64), ${k1}$) + mean = sum(x, mask) / n + + #:if t1[0] == 'r' + res = sum((x - mean)**2, mask) / (n - 1._${k1}$) + #:else + res = sum(abs(x - mean)**2, mask) / (n - 1._${k1}$) + #:endif + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask_all",rank, t1, k1, 'dp') + module function ${RName}$(x, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res + + real(dp) :: n, mean + + n = real(count(mask, kind = int64), dp) + mean = sum(real(x, dp), mask) / n + + res = sum((real(x, dp) - mean)**2, mask) / (n - 1._dp) + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in RC_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask",rank, t1, k1) + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(${k1}$) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$ + ${t1}$ :: mean${reduced_shape('x', rank, 'dim')}$ + + res = 0._${k1}$ + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(count(mask, dim), ${k1}$) + mean = sum(x, dim, mask) / n + do i = 1, size(x, dim) + #:if t1[0] == 'r' + res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2,& + #:else + res = res + merge( abs(x${select_subarray(rank, [(fi, 'i')])}$ - mean)**2,& + #:endif + 0._${k1}$,& + mask${select_subarray(rank, [(fi, 'i')])}$) + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._${k1}$) + + end function ${RName}$ + #:endfor + #:endfor + + + #:for k1, t1 in INT_KINDS_TYPES + #:for rank in RANKS + #:set RName = rname("var_mask",rank, t1, k1, 'dp') + module function ${RName}$(x, dim, mask) result(res) + ${t1}$, intent(in) :: x${ranksuffix(rank)}$ + integer, intent(in) :: dim + logical, intent(in) :: mask${ranksuffix(rank)}$ + real(dp) :: res${reduced_shape('x', rank, 'dim')}$ + + integer :: i + real(dp) :: n${reduced_shape('x', rank, 'dim')}$ + real(dp) :: mean${reduced_shape('x', rank, 'dim')}$ + + res = 0._dp + select case(dim) + #:for fi in range(1, rank+1) + case(${fi}$) + n = real(count(mask, dim), dp) + mean = sum(real(x, dp), dim, mask) / n + do i = 1, size(x, dim) + res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean)**2,& + 0._dp, mask${select_subarray(rank, [(fi, 'i')])}$) + end do + #:endfor + case default + call error_stop("ERROR (mean): wrong dimension") + end select + res = res / (n - 1._dp) + + end function ${RName}$ + #:endfor + #:endfor + +end submodule diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index e86fe23d7..239ac2373 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -1,4 +1,5 @@ ADDTEST(mean) +ADDTEST(var) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 9faf154cb..1d25a5509 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,3 +1,3 @@ -PROGS_SRC = test_mean.f90 +PROGS_SRC = test_mean.f90 test_var.f90 include ../Makefile.manual.test.mk diff --git a/src/tests/stats/test_var.f90 b/src/tests/stats/test_var.f90 new file mode 100644 index 000000000..868c7a194 --- /dev/null +++ b/src/tests/stats/test_var.f90 @@ -0,0 +1,478 @@ +program test_var + use stdlib_experimental_error, only: assert + use stdlib_experimental_kinds, only: sp, dp, int32, int64 + use stdlib_experimental_io, only: loadtxt + use stdlib_experimental_stats, only: mean, var + implicit none + + + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) + real(dp), parameter :: dptol = 1000 * epsilon(1._dp) + + integer(int32) :: i321(5) = [1, 2, 3, 4, 5] + integer(int64) :: i641(5) = [1, 2, 3, 4, 5] + + integer(int32), allocatable :: i32(:,:), i323(:, :, :) + integer(int64), allocatable :: i64(:,:), i643(:, :, :) + + real(sp) :: s1(5) = [1.0_sp, 2.0_sp, 3.0_sp, 4.0_sp, 5.0_sp] + real(dp) :: d1(5) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp] + + real(sp), allocatable :: s(:, :), s3(:, :, :) + real(dp), allocatable :: d3(:, :, :) + real(dp) :: d(4, 3) = reshape([1._dp, 3._dp, 5._dp, 7._dp,& + 2._dp, 4._dp, 6._dp, 8._dp,& + 9._dp, 10._dp, 11._dp, 12._dp], [4, 3]) + + + complex(sp) :: cs1(5) = [ cmplx(0.57706_sp, 0.00000_sp),& + cmplx(0.00000_sp, 1.44065_sp),& + cmplx(1.26401_sp, 0.00000_sp),& + cmplx(0.00000_sp, 0.88833_sp),& + cmplx(1.14352_sp, 0.00000_sp)] + complex(dp) :: cd1(5) = [ cmplx(0.57706_dp, 0.00000_dp),& + cmplx(0.00000_dp, 1.44065_dp),& + cmplx(1.26401_dp, 0.00000_dp),& + cmplx(0.00000_dp, 0.88833_dp),& + cmplx(1.14352_dp, 0.00000_dp)] + complex(sp) :: cs(5,3) + complex(dp) :: cd(5,3) + + + !sp + !1dim + print*,' test_sp_1dim' + call assert( abs(var(s1) - 2.5) < sptol) + call assert( abs(var(s1, dim=1) - 2.5) < sptol) + + print*,' test_sp_1dim_mask' + call assert( isnan(var(s1, .false.))) + call assert( isnan(var(s1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(s1, s1 < 5) - 5./3.) < sptol) + call assert( abs(var(s1, 1, s1 < 5) - 5./3.) < sptol) + + + !2dim + print*,' test_sp_2dim' + s = d + call assert( abs(var(s) - 13) < sptol) + call assert( all( abs( var(s, 1) - [20. / 3., 20. / 3., 5. / 3.]) < sptol)) + call assert( all( abs( var(s, 2) - [19.0, 43. / 3., 31. / 3. , 7.0]) < sptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(s, .false.))) + call assert( any(isnan(var(s, 1, .false.)))) + call assert( any(isnan(var(s, 2, .false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(s, s < 11) - 27.5 / 3.) < sptol) + call assert( all( abs( var(s, 1, s < 11) - [20. / 3., 20. / 3., 0.5]) < sptol)) + call assert( all( abs( var(s, 2, s < 11) - [19.0, 43. / 3., 0.5 , 0.5]) < sptol)) + + + !3dim + allocate(s3(size(s,1),size(s,2),3)) + s3(:,:,1)=s; + s3(:,:,2)=s*2; + s3(:,:,3)=s*4; + + print*,' test_sp_3dim' + call assert( abs(var(s3) - 153.4) < sptol) + call assert( all( abs( var(s3, 1) -& + reshape([20. / 3., 20. / 3., 5. / 3.,& + 4* 20. / 3., 4* 20. / 3., 4* 5. / 3.,& + 16* 20. / 3., 16* 20. / 3., 16* 5. / 3.],& + [size(s3,2), size(s3,3)]))& + < sptol)) + call assert( all( abs( var(s3, 2) -& + reshape([19.0, 43. / 3., 31. / 3. , 7.0,& + 4* 19.0, 4* 43. / 3., 4* 31. / 3. , 4* 7.0,& + 16* 19.0, 16* 43. / 3., 16* 31. / 3. , 16* 7.0],& + [size(s3,1), size(s3,3)] ))& + < sptol)) + call assert( all(abs( var(s3, 3) -& + reshape([ 7./3., 21., 175./3.,& + 343./3., 28./3., 112./3.,& + 84., 448./3., 189.,& + 700./3., 847./3., 336.], [size(s3,1), size(s3,2)] ))& + < sptol)) + + print*,' test_sp_3dim_mask' + call assert( isnan(var(s3, .false.))) + call assert( any(isnan(var(s3, 1, .false.)))) + call assert( any(isnan(var(s3, 2, .false.)))) + call assert( any(isnan(var(s3, 3, .false.)))) + + print*,' test_sp_3dim_mask_array' + call assert( abs(var(s3, s3 < 11) - 8.2205877_sp) < sptol) + call assert( all( abs( var(s3, 1, s3 < 25) -& + reshape([20./3., 20./3., 5./3., 80./3., 80./3., 20./3., 64., 64., 0.],& + [size(s3, 2), size(s3, 3)])) < sptol )) + call assert( all( abs( var(s3, 2, s3 < 25) -& + reshape([19., 43./3., 31./3., 7.,& + 4*19., 4*43./3., 4*31./3., 4*7.,& + 8., 8., 8., 0.],& + [size(s3, 1), size(s3, 3)])) < sptol )) + call assert( all( abs( var(s3, 3, s3 < 25) -& + reshape([ 7./3., 21., 175./3.,& + 24.5, 28./3., 112./3.,& + 84., 32., 40.5,& + 50., 60.5, 72.], [size(s3,1), size(s3,2)] ))& + < sptol )) + + + !dp + !1dim + print*,' test_dp_1dim' + call assert( abs(var(d1) - 2.5) < dptol) + call assert( abs(var(d1, 1) - 2.5) < dptol) + + print*,' test_dp_1dim_mask' + call assert( isnan(var(d1, .false.))) + call assert( isnan(var(d1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(d1, d1 < 5) - 5._dp/3._dp) < dptol) + call assert( abs(var(d1, 1, d1 < 5) - 5._dp/3._dp) < dptol) + + !2dim + print*,' test_dp_2dim' + call assert( abs(var(d) - 13) < dptol) + call assert( all( abs( var(d,1) -& + [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol)) + call assert( all( abs( var(d,2) -& + [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol)) + + print*,' test_dp_2dim_mask' + call assert( isnan(var(d, .false.))) + call assert( any(isnan(var(d, 1, .false.)))) + call assert( any(isnan(var(d, 2, .false.)))) + + print*,' test_dp_2dim_mask_array' + call assert( abs(var(d, d < 11) - 27.5_dp / 3._dp) < dptol) + call assert( all( abs( var(d, 1, d < 11) -& + [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol)) + call assert( all( abs( var(d, 2, d < 11) -& + [19.0_dp, 43._dp / 3._dp, 0.5 _dp, 0.5_dp]) < dptol)) + + !3dim + allocate(d3(size(d,1),size(d,2),3)) + d3(:,:,1)=d; + d3(:,:,2)=d*2; + d3(:,:,3)=d*4; + + print*,' test_dp_3dim' + call assert( abs(var(d3) - 153.4_dp) < dptol) + call assert( all( abs( var(d3, 1) -& + reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,& + 4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,& + 16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],& + [size(d3,2), size(d3,3)]))& + < dptol)) + print*,' test_dp_3dim' + call assert( all( abs( var(d3, 2) -& + reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,& + 16* 7.0_dp],& + [size(d3,1), size(d3,3)] ))& + < dptol)) + print*,' test_dp_3dim' + call assert( all(abs( var(d3, 3) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 448._dp/3._dp, 189._dp,& + 700._dp/3._dp, 847._dp/3._dp, 336._dp],& + [size(d3,1), size(d3,2)] ))& + < dptol)) + + print*,' test_dp_3dim_mask' + call assert( isnan(var(d3, .false.))) + call assert( any(isnan(var(d3, 1, .false.)))) + call assert( any(isnan(var(d3, 2, .false.)))) + call assert( any(isnan(var(d3, 3, .false.)))) + + print*,' test_dp_3dim_mask_array' + call assert( abs(var(d3, d3 < 25) - 46.041379310344823_dp) < dptol) + call assert( all( abs( var(d3, 1, d3 < 25) -& + reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,& + 80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,& + 64._dp, 64._dp, 0._dp],& + [size(d3, 2), size(d3, 3)]))& + < dptol )) + call assert( all( abs( var(d3, 2, d3 < 25) -& + reshape([19._dp, 43._dp/3._dp, 31._dp/3._dp, 7._dp,& + 4*19._dp, 4*43._dp/3._dp, 4*31._dp/3._dp, 4*7._dp,& + 8._dp, 8._dp, 8._dp, 0._dp],& + [size(d3, 1), size(d3, 3)]))& + < dptol )) + call assert( all( abs( var(d3, 3, d3 < 25) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 24.5_dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 32._dp, 40.5_dp,& + 50._dp, 60.5_dp, 72._dp],& + [size(d3,1), size(d3,2)] ))& + < dptol )) + + + + !int32 + !1dim + print*,' test_int32_1dim' + call assert( abs(var(i321) - 2.5) < dptol) + call assert( abs(var(i321, 1) - 2.5) < dptol) + + print*,' test_int32_1dim_mask' + call assert( isnan(var(i321, .false.))) + call assert( isnan(var(i321, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(i321, i321 < 5) - 5._dp/3._dp) < dptol) + call assert( abs(var(i321, 1, i321 < 5) - 5._dp/3._dp) < dptol) + + !2dim + print*,' test_int32_2dim' + i32 = d + call assert( abs(var(i32) - 13) < dptol) + call assert( all( abs( var(i32,1) -& + [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol)) + call assert( all( abs( var(i32,2) -& + [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol)) + + print*,' test_int32_2dim_mask' + call assert( isnan(var(i32, .false.))) + call assert( any(isnan(var(i32, 1, .false.)))) + call assert( any(isnan(var(i32, 2, .false.)))) + + print*,' test_int32_2dim_mask_array' + call assert( abs(var(i32, i32 < 11) - 27.5_dp / 3._dp) < dptol) + call assert( all( abs( var(i32, 1, i32 < 11) -& + [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol)) + call assert( all( abs( var(i32, 2, i32 < 11) -& + [19.0_dp, 43._dp / 3._dp, 0.5_dp, 0.5_dp]) < dptol)) + + !3dim + allocate(i323(size(d,1),size(d,2),3)) + i323(:,:,1)=d; + i323(:,:,2)=d*2; + i323(:,:,3)=d*4; + + print*,' test_int32_3dim' + call assert( abs(var(i323) - 153.4_dp) < dptol) + call assert( all( abs( var(i323, 1) -& + reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,& + 4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,& + 16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],& + [size(i323,2), size(i323,3)]))& + < dptol)) + call assert( all( abs( var(i323, 2) -& + reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,& + 16* 7.0_dp],& + [size(i323,1), size(i323,3)] ))& + < dptol)) + call assert( all(abs( var(i323, 3) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 448._dp/3._dp, 189._dp,& + 700._dp/3._dp, 847._dp/3._dp, 336._dp],& + [size(i323,1), size(i323,2)] ))& + < dptol)) + + print*,' test_int32_3dim_mask' + call assert( isnan(var(i323, .false.))) + call assert( any(isnan(var(i323, 1, .false.)))) + call assert( any(isnan(var(i323, 2, .false.)))) + call assert( any(isnan(var(i323, 3, .false.)))) + + print*,' test_int32_3dim_mask_array' + call assert( abs(var(i323, i323 < 25) - 46.041379310344823_dp) < dptol) + call assert( all( abs( var(i323, 1, i323 < 25) -& + reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,& + 80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,& + 64._dp, 64._dp, 0._dp],& + [size(i323, 2), size(i323, 3)]))& + < dptol )) + call assert( all( abs( var(i323, 2, i323 < 25) -& + reshape([19._dp, 43._dp/3._dp, 31._dp/3._dp, 7._dp,& + 4*19._dp, 4*43._dp/3._dp, 4*31._dp/3._dp, 4*7._dp,& + 8._dp, 8._dp, 8._dp, 0._dp],& + [size(i323, 1), size(i323, 3)]))& + < dptol )) + call assert( all( abs( var(i323, 3, i323 < 25) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 24.5_dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 32._dp, 40.5_dp,& + 50._dp, 60.5_dp, 72._dp],& + [size(i323,1), size(i323,2)] ))& + < dptol )) + + + !int64 + !1dim + print*,' test_int64_1dim' + call assert( abs(var(i641) - 2.5) < dptol) + call assert( abs(var(i641, 1) - 2.5) < dptol) + + print*,' test_int64_1dim_mask' + call assert( isnan(var(i641, .false.))) + call assert( isnan(var(i641, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(i641, i641 < 5) - 5._dp/3._dp) < dptol) + call assert( abs(var(i641, 1, i641 < 5) - 5._dp/3._dp) < dptol) + + !2dim + print*,' test_int64_2dim' + i64 = d + call assert( abs(var(i64) - 13) < dptol) + call assert( all( abs( var(i64,1) -& + [20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp]) < dptol)) + call assert( all( abs( var(i64,2) -& + [19.0_dp, 43._dp/3._dp, 31._dp/3._dp, 7.0_dp]) < dptol)) + + print*,' test_int64_2dim_mask' + call assert( isnan(var(i64, .false.))) + call assert( any(isnan(var(i64, 1, .false.)))) + call assert( any(isnan(var(i64, 2, .false.)))) + + print*,' test_int64_2dim_mask_array' + call assert( abs(var(i64, i64 < 11) - 27.5_dp / 3._dp) < dptol) + call assert( all( abs( var(i64, 1, i64 < 11) -& + [20._dp / 3._dp, 20._dp / 3._dp, 0.5_dp]) < dptol)) + call assert( all( abs( var(i64, 2, i64 < 11) -& + [19.0_dp, 43._dp / 3._dp, 0.5 _dp, 0.5_dp]) < dptol)) + + !3dim + allocate(i643(size(d,1),size(d,2),3)) + i643(:,:,1)=d; + i643(:,:,2)=d*2; + i643(:,:,3)=d*4; + + print*,' test_int32_3dim' + call assert( abs(var(i643) - 153.4_dp) < dptol) + call assert( all( abs( var(i643, 1) -& + reshape([20._dp / 3._dp, 20._dp / 3._dp, 5._dp / 3._dp,& + 4* 20._dp / 3._dp, 4* 20._dp / 3._dp, 4* 5._dp / 3._dp,& + 16* 20._dp / 3._dp, 16* 20._dp / 3._dp, 16* 5._dp / 3._dp],& + [size(i643,2), size(i643,3)]))& + < dptol)) + call assert( all( abs( var(i643, 2) -& + reshape([19.0_dp, 43._dp / 3._dp, 31._dp / 3._dp , 7.0_dp,& + 4* 19.0_dp, 4* 43._dp / 3._dp, 4* 31._dp / 3._dp , 4* 7.0_dp,& + 16* 19.0_dp, 16* 43._dp / 3._dp, 16* 31._dp / 3._dp ,& + 16* 7.0_dp],& + [size(i643,1), size(i643,3)] ))& + < dptol)) + call assert( all(abs( var(i643, 3) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 343._dp/3._dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 448._dp/3._dp, 189._dp,& + 700._dp/3._dp, 847._dp/3._dp, 336._dp],& + [size(i643,1), size(i643,2)] ))& + < dptol)) + + print*,' test_int32_3dim_mask' + call assert( isnan(var(i643, .false.))) + call assert( any(isnan(var(i643, 1, .false.)))) + call assert( any(isnan(var(i643, 2, .false.)))) + call assert( any(isnan(var(i643, 3, .false.)))) + + print*,' test_int32_3dim_mask_array' + call assert( abs(var(i643, i643 < 25) - 46.041379310344823_dp) < dptol) + call assert( all( abs( var(i643, 1, i643 < 25) -& + reshape([20._dp/3._dp, 20._dp/3._dp, 5._dp/3._dp,& + 80._dp/3._dp, 80._dp/3._dp, 20._dp/3._dp,& + 64._dp, 64._dp, 0._dp],& + [size(i643, 2), size(i643, 3)]))& + < dptol )) + call assert( all( abs( var(i643, 2, i643 < 25) -& + reshape([19._dp, 43._dp/3._dp, 31._dp/3._dp, 7._dp,& + 4*19._dp, 4*43._dp/3._dp, 4*31._dp/3._dp, 4*7._dp,& + 8._dp, 8._dp, 8._dp, 0._dp],& + [size(i643, 1), size(i643, 3)]))& + < dptol )) + call assert( all( abs( var(i643, 3, i643 < 25) -& + reshape([ 7._dp/3._dp, 21._dp, 175._dp/3._dp,& + 24.5_dp, 28._dp/3._dp, 112._dp/3._dp,& + 84._dp, 32._dp, 40.5_dp,& + 50._dp, 60.5_dp, 72._dp],& + [size(i643,1), size(i643,2)] ))& + < dptol )) + + !csp + !1dim + print*,' test_csp_1dim' + call assert( abs(var(cs1) - (var(real(cs1)) + var(aimag(cs1)))) < sptol) + call assert( abs(var(cs1, dim=1) - (var(real(cs1),1) + var(aimag(cs1), 1)) ) < sptol) + + print*,' test_csp_1dim_mask' + call assert( isnan(var(cs1, .false.))) + call assert( isnan(var(cs1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(cs1, aimag(cs1) == 0) - var(real(cs1), aimag(cs1) == 0)) < sptol) + call assert( abs(var(cs1, 1, aimag(cs1) == 0) - var(real(cs1), 1, aimag(cs1) == 0)) < sptol) + + !2dim + cs(:,1) = cs1 + cs(:,2) = cs1*3_sp + cs(:,3) = cs1*1.5_sp + + print*,' test_sp_2dim' + print*,var(cs),(var(real(cs)) + var(aimag(cs))) + print*,var(cs)-(var(real(cs)) + var(aimag(cs))),sptol + call assert( abs(var(cs) - (var(real(cs)) + var(aimag(cs)))) < sptol) + call assert( all( abs( var(cs, 1) - (var(real(cs), 1) + var(aimag(cs), 1))) < sptol)) + call assert( all( abs( var(cs, 2) - (var(real(cs), 2) + var(aimag(cs), 2))) < sptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(cs, .false.))) + call assert( any(isnan(var(cs, 1, .false.)))) + call assert( any(isnan(var(cs, 2, .false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(cs, aimag(cs) == 0) - var(real(cs), aimag(cs) == 0)) < sptol) + call assert( all( abs( var(cs, 1, aimag(cs) == 0) - var(real(cs), 1, aimag(cs) == 0)) < sptol)) + call assert( all( abs( var(cs, 2, aimag(cs) == 0) - var(real(cs), 2, aimag(cs) == 0)) < sptol)) + + !cdp + !1dim + print*,' test_cdp_1dim' + call assert( abs(var(cd1) - (var(real(cd1)) + var(aimag(cd1)))) < dptol) + call assert( abs(var(cd1, dim=1) - (var(real(cd1),1) + var(aimag(cd1), 1)) ) < dptol) + + print*,' test_cdp_1dim_mask' + call assert( isnan(var(cd1, .false.))) + call assert( isnan(var(cd1, 1, .false.))) + + print*,' test_sp_1dim_mask_array' + call assert( abs(var(cd1, aimag(cd1) == 0) - var(real(cd1), aimag(cd1) == 0)) < dptol) + call assert( abs(var(cd1, 1, aimag(cd1) == 0) - var(real(cd1), 1, aimag(cd1) == 0)) < dptol) + + !2dim + cd(:,1) = cd1 + cd(:,2) = cd1*3_sp + cd(:,3) = cd1*1.5_sp + + print*,' test_sp_2dim' + print*,var(cd),(var(real(cd)) + var(aimag(cd))) + print*,var(cd)-(var(real(cd)) + var(aimag(cd))),dptol + call assert( abs(var(cd) - (var(real(cd)) + var(aimag(cd)))) < dptol) + call assert( all( abs( var(cd, 1) - (var(real(cd), 1) + var(aimag(cd), 1))) < dptol)) + call assert( all( abs( var(cd, 2) - (var(real(cd), 2) + var(aimag(cd), 2))) < dptol)) + + print*,' test_sp_2dim_mask' + call assert( isnan(var(cd, .false.))) + call assert( any(isnan(var(cd, 1, .false.)))) + call assert( any(isnan(var(cd, 2, .false.)))) + + print*,' test_sp_2dim_mask_array' + call assert( abs(var(cd, aimag(cd) == 0) - var(real(cd), aimag(cd) == 0)) < dptol) + call assert( all( abs( var(cd, 1, aimag(cd) == 0) - var(real(cd), 1, aimag(cd) == 0)) < dptol)) + call assert( all( abs( var(cd, 2, aimag(cd) == 0) - var(real(cd), 2, aimag(cd) == 0)) < dptol)) + +end program