From ae868d3ada099a6f1a15ed711b72d2ffa2637df9 Mon Sep 17 00:00:00 2001 From: Zhuowen Zhao Date: Thu, 14 Sep 2017 15:25:22 -0400 Subject: [PATCH 1/2] only one error for math_check; new "math_expand('what' by 'how')" --- src/IO.f90 | 8 +---- src/math.f90 | 95 ++++++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 79 insertions(+), 24 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 3d80feefa..d067a84c0 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -1557,13 +1557,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (400_pInt) msg = 'matrix inversion error' case (401_pInt) - msg = 'math_check: quat -> axisAngle -> quat failed' - case (402_pInt) - msg = 'math_check: quat -> R -> quat failed' - case (403_pInt) - msg = 'math_check: quat -> euler -> quat failed' - case (404_pInt) - msg = 'math_check: R -> euler -> R failed' + msg = 'math_check failed' case (405_pInt) msg = 'I_TO_HALTON-error: an input base BASE is <= 1' case (406_pInt) diff --git a/src/math.f90 b/src/math.f90 index b5507d868..1bf903ced 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -180,15 +180,11 @@ subroutine math_init implicit none integer(pInt) :: i - real(pReal), dimension(3,3) :: R,R2 - real(pReal), dimension(3) :: Eulers,v - real(pReal), dimension(4) :: q,q2,axisangle,randTest + real(pReal), dimension(4) :: randTest ! the following variables are system dependend and shound NOT be pInt integer :: randSize ! gfortran requires a variable length to compile integer, dimension(:), allocatable :: randInit ! if recalculations of former randomness (with given seed) is necessary ! comment the first random_seed call out, set randSize to 1, and use ifort - character(len=64) :: error_msg - write(6,'(/,a)') ' <<<+- math init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -211,7 +207,7 @@ subroutine math_init enddo write(6,'(a,I2)') ' size of random seed: ', randSize - do i =1, randSize + do i = 1_pInt,randSize write(6,'(a,I2,I14)') ' value of random seed: ', i, randInit(i) enddo write(6,'(a,4(/,26x,f17.14),/)') ' start of random sequence: ', randTest @@ -221,6 +217,27 @@ subroutine math_init call halton_seed_set(int(randInit(1), pInt)) call halton_ndim_set(3_pInt) + call math_check() + +end subroutine math_init + +!-------------------------------------------------------------------------------------------------- +!> @brief check correctness of (some) math functions +!-------------------------------------------------------------------------------------------------- +subroutine math_check + + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use prec, only: tol_math_check + use numerics, only: & + fixedSeed + use IO, only: IO_error, IO_timeStamp + + character(len=64) :: error_msg + + real(pReal), dimension(3,3) :: R,R2 + real(pReal), dimension(3) :: Eulers,v + real(pReal), dimension(4) :: q,q2,axisangle + ! --- check rotation dictionary --- q = math_qRand() ! random quaternion @@ -230,7 +247,8 @@ subroutine math_init q2 = math_axisAngleToQ(axisangle(1:3),axisangle(4)) if ( any(abs( q-q2) > tol_math_check) .and. & any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) + write (error_msg, '(a,e14.6)' ) & + 'quat -> axisAngle -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) call IO_error(401_pInt,ext_msg=error_msg) endif @@ -239,8 +257,9 @@ subroutine math_init q2 = math_RtoQ(R) if ( any(abs( q-q2) > tol_math_check) .and. & any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(402_pInt,ext_msg=error_msg) + write (error_msg, '(a,e14.6)' ) & + 'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) + call IO_error(401_pInt,ext_msg=error_msg) endif ! +++ q -> euler -> q +++ @@ -248,28 +267,46 @@ subroutine math_init q2 = math_EulerToQ(Eulers) if ( any(abs( q-q2) > tol_math_check) .and. & any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) 'maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(403_pInt,ext_msg=error_msg) + write (error_msg, '(a,e14.6)' ) & + 'quat -> euler -> quatmaximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) + call IO_error(401_pInt,ext_msg=error_msg) endif ! +++ R -> euler -> R +++ Eulers = math_RtoEuler(R) R2 = math_EulerToR(Eulers) if ( any(abs( R-R2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) 'maximum deviation ',maxval(abs( R-R2)) - call IO_error(404_pInt,ext_msg=error_msg) + write (error_msg, '(a,e14.6)' ) & + 'R -> euler -> R maximum deviation ',maxval(abs( R-R2)) + call IO_error(401_pInt,ext_msg=error_msg) endif ! +++ check rotation sense of q and R +++ - q = math_qRand() ! random quaternion call halton(3_pInt,v) ! random vector R = math_qToR(q) if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then - write(6,'(a,4(f8.3,1x))') 'q',q - call IO_error(409_pInt) + write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v' + call IO_error(401_pInt,ext_msg=error_msg) endif -end subroutine math_init + ! +++ check vector expansion +++ + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & + math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt,3_pInt,0_pInt])) > tol_math_check)) then + write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]' + call IO_error(401_pInt,ext_msg=error_msg) + endif + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal] - & + math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1_pInt,2_pInt])) > tol_math_check)) then + write (error_msg, '(a)' ) 'math_expand [1,2,3] by [1,2] => [1,2,2]' + call IO_error(401_pInt,ext_msg=error_msg) + endif + if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,1.0_pReal,1.0_pReal,1.0_pReal] - & + math_expand([1.0_pReal,2.0_pReal],[1_pInt,2_pInt,3_pInt])) > tol_math_check)) then + write (error_msg, '(a)' ) 'math_expand [1,2] by [1,2,3] => [1,2,2,1,1,1]' + call IO_error(401_pInt,ext_msg=error_msg) + endif + +end subroutine math_check !-------------------------------------------------------------------------------------------------- @@ -334,6 +371,30 @@ recursive subroutine math_qsort(a, istart, iend) end subroutine math_qsort +!-------------------------------------------------------------------------------------------------- +!> @brief vector expansion +!> @details takes a set of numbers (a,b,c,...) and corresponding multipliers (x,y,z,...) +!> to return a vector of x times a, y times b, z times c, ... +!-------------------------------------------------------------------------------------------------- +pure function math_expand(what,how) + + implicit none + real(pReal), dimension(:), intent(in) :: what + integer(pInt), dimension(:), intent(in) :: how + real(pReal), dimension(sum(how)) :: math_expand + integer(pInt) :: i,j,o + + o = 0_pInt + do i = 1, size(how) + do j = 1, how(i) + o = o + 1_pInt + math_expand(o) = what(1+mod(i-1,size(what))) + enddo + enddo + +end function math_expand + + !-------------------------------------------------------------------------------------------------- !> @brief range of integers starting at one !-------------------------------------------------------------------------------------------------- From fe4952b100085b730fad5ccb012e5072344106a1 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 15 Sep 2017 06:20:35 +0200 Subject: [PATCH 2/2] [skip ci] updated version information after successful test of v2.0.1-918-gae868d3 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 229e649b1..3cad8e739 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.1-909-g9e502c7 +v2.0.1-918-gae868d3