only one error for math_check; new "math_expand('what' by 'how')"
This commit is contained in:
Zhuowen Zhao
parent
e211fa614a
commit
ae868d3ada
|
|
@ -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)
|
||||
|
|
|
|||
95
src/math.f90
95
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
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
Loading…
Reference in New Issue