2012-08-09 16:31:53 +05:30
!--------------------------------------------------------------------------------------------------
2013-09-14 16:29:35 +05:30
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2014-05-08 20:25:19 +05:30
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
2016-03-21 03:50:58 +05:30
!> @brief setting precision for real and int type
2012-08-09 16:31:53 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-06 20:22:48 +05:30
module prec
2019-05-11 01:18:43 +05:30
use , intrinsic :: IEEE_arithmetic
2020-09-19 14:20:32 +05:30
use , intrinsic :: ISO_C_binding
2019-03-07 23:59:48 +05:30
2019-03-06 19:54:42 +05:30
implicit none
2019-09-23 12:23:56 +05:30
public
2020-02-29 19:30:47 +05:30
2019-03-06 19:54:42 +05:30
! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
2020-02-29 19:30:47 +05:30
integer , parameter :: pReal = IEEE_selected_real_kind ( 15 , 307 ) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
2020-09-12 12:44:33 +05:30
integer , parameter :: pI32 = selected_int_kind ( 9 ) !< number with at least up to +-1e9 (typically 32 bit)
integer , parameter :: pI64 = selected_int_kind ( 18 ) !< number with at least up to +-1e18 (typically 64 bit)
2019-03-06 20:22:52 +05:30
#if(INT==8)
2020-09-12 12:44:33 +05:30
integer , parameter :: pInt = pI64
2012-08-31 01:56:28 +05:30
#else
2020-09-12 12:44:33 +05:30
integer , parameter :: pInt = pI32
2012-08-28 21:38:17 +05:30
#endif
2020-02-29 19:30:47 +05:30
integer , parameter :: pStringLen = 256 !< default string length
integer , parameter :: pPathLen = 4096 !< maximum length of a path name on linux
2019-03-06 19:54:42 +05:30
2020-02-29 19:30:47 +05:30
real ( pReal ) , parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
2019-03-06 19:54:42 +05:30
2020-02-29 19:30:47 +05:30
type :: group_float !< variable length datatype used for storage of state
2019-03-06 19:54:42 +05:30
real ( pReal ) , dimension ( : ) , pointer :: p
end type group_float
! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
2020-02-29 19:30:47 +05:30
type :: tState
2019-03-07 23:59:48 +05:30
integer :: &
sizeState = 0 , & !< size of state
sizeDotState = 0 , & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
offsetDeltaState = 0 , & !< index offset of delta state
2019-12-11 04:40:02 +05:30
sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
2019-03-06 19:54:42 +05:30
real ( pReal ) , pointer , dimension ( : ) , contiguous :: &
2020-03-15 14:21:40 +05:30
atol
2019-03-07 23:59:48 +05:30
real ( pReal ) , pointer , dimension ( : , : ) , contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous
2019-03-06 19:54:42 +05:30
state0 , &
state , & !< state
dotState , & !< rate of state change
deltaState !< increment of state change
real ( pReal ) , allocatable , dimension ( : , : ) :: &
partionedState0 , &
2020-04-01 01:32:53 +05:30
subState0
2019-03-06 19:54:42 +05:30
end type
2020-02-29 19:30:47 +05:30
type , extends ( tState ) :: tPlasticState
2019-03-07 23:59:48 +05:30
logical :: &
2019-03-06 19:54:42 +05:30
nonlocal = . false .
real ( pReal ) , pointer , dimension ( : , : ) :: &
2020-04-01 13:30:12 +05:30
slipRate !< slip rate
2019-03-06 19:54:42 +05:30
end type
2020-02-29 19:30:47 +05:30
type :: tSourceState
2019-03-06 19:54:42 +05:30
type ( tState ) , dimension ( : ) , allocatable :: p !< tState for each active source mechanism in a phase
end type
2019-03-07 23:59:48 +05:30
2020-02-29 19:30:47 +05:30
type :: tHomogMapping
2019-03-07 23:59:48 +05:30
integer , pointer , dimension ( : , : ) :: p
end type
2019-03-06 19:54:42 +05:30
2019-03-07 23:59:48 +05:30
real ( pReal ) , private , parameter :: PREAL_EPSILON = epsilon ( 0.0_pReal ) !< minimum positive number such that 1.0 + EPSILON /= 1.0.
2019-03-07 15:32:27 +05:30
real ( pReal ) , private , parameter :: PREAL_MIN = tiny ( 0.0_pReal ) !< smallest normalized floating point number
2019-03-06 19:54:42 +05:30
2020-02-29 19:30:47 +05:30
integer , dimension ( 0 ) , parameter :: &
2020-01-02 23:14:51 +05:30
emptyIntArray = [ integer :: ]
2020-02-29 19:30:47 +05:30
real ( pReal ) , dimension ( 0 ) , parameter :: &
2020-01-02 23:14:51 +05:30
emptyRealArray = [ real ( pReal ) :: ]
2020-02-29 19:30:47 +05:30
character ( len = pStringLen ) , dimension ( 0 ) , parameter :: &
2020-01-02 23:14:51 +05:30
emptyStringArray = [ character ( len = pStringLen ) :: ]
2019-12-21 15:49:33 +05:30
2019-09-23 12:23:56 +05:30
private :: &
2020-05-16 20:35:03 +05:30
selfTest
2019-03-07 23:59:48 +05:30
2012-03-06 20:22:48 +05:30
contains
2012-08-28 21:38:17 +05:30
2015-04-13 15:32:52 +05:30
2012-08-09 16:31:53 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-06 21:06:05 +05:30
!> @brief report precision and do self test
2012-08-09 16:31:53 +05:30
!--------------------------------------------------------------------------------------------------
2012-03-06 20:22:48 +05:30
subroutine prec_init
2015-04-13 15:32:52 +05:30
2020-09-17 22:58:41 +05:30
print '(/,a)' , ' <<<+- prec init -+>>>'
print '(a,i3)' , ' Size of integer in bit: ' , bit_size ( 0 )
print '(a,i19)' , ' Maximum value: ' , huge ( 0 )
print '(/,a,i3)' , ' Size of float in bit: ' , storage_size ( 0.0_pReal )
print '(a,e10.3)' , ' Maximum value: ' , huge ( 0.0_pReal )
print '(a,e10.3)' , ' Minimum value: ' , tiny ( 0.0_pReal )
print '(a,i3)' , ' Decimal precision: ' , precision ( 0.0_pReal )
2014-04-15 14:50:38 +05:30
2020-05-16 20:35:03 +05:30
call selfTest
2019-03-07 23:59:48 +05:30
2012-03-06 20:22:48 +05:30
end subroutine prec_init
2009-08-31 20:39:15 +05:30
2015-04-13 15:32:52 +05:30
2016-03-06 02:07:22 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Test floating point numbers with double precision for equality.
2016-03-21 03:50:58 +05:30
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
2016-10-18 23:20:05 +05:30
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
2017-11-21 13:48:03 +05:30
! AlmostEqualRelative
2016-03-06 02:07:22 +05:30
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq ( a , b , tol )
2016-05-25 11:22:56 +05:30
2019-03-06 20:27:42 +05:30
real ( pReal ) , intent ( in ) :: a , b
real ( pReal ) , intent ( in ) , optional :: tol
real ( pReal ) :: eps
2019-03-07 23:59:48 +05:30
2019-03-06 20:27:42 +05:30
if ( present ( tol ) ) then
eps = tol
else
2019-03-07 15:32:27 +05:30
eps = PREAL_EPSILON * maxval ( abs ( [ a , b ] ) )
2019-03-06 20:27:42 +05:30
endif
2019-09-23 21:59:38 +05:30
dEq = merge ( . True . , . False . , abs ( a - b ) < = eps )
2016-05-25 11:22:56 +05:30
2016-03-06 02:07:22 +05:30
end function dEq
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Test floating point numbers with double precision for inequality.
2016-03-21 03:50:58 +05:30
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
2016-10-18 23:20:05 +05:30
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
2017-11-21 13:48:03 +05:30
! AlmostEqualRelative NOT
2016-03-06 02:07:22 +05:30
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq ( a , b , tol )
2016-05-25 11:22:56 +05:30
2019-03-06 20:27:42 +05:30
real ( pReal ) , intent ( in ) :: a , b
real ( pReal ) , intent ( in ) , optional :: tol
2019-03-07 23:59:48 +05:30
2019-03-06 20:27:42 +05:30
if ( present ( tol ) ) then
2019-09-23 12:23:56 +05:30
dNeq = . not . dEq ( a , b , tol )
2019-03-06 20:27:42 +05:30
else
2019-09-23 12:23:56 +05:30
dNeq = . not . dEq ( a , b )
2019-03-06 20:27:42 +05:30
endif
2016-03-06 02:07:22 +05:30
end function dNeq
2016-06-30 14:00:40 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Test floating point number with double precision for equality to 0.
2017-11-21 13:48:03 +05:30
! replaces "==0" but everything not representable as a normal number is treated as 0. Counterpart to dNeq0
! https://de.mathworks.com/help/matlab/ref/realmin.html
! https://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
2016-06-30 14:00:40 +05:30
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq0 ( a , tol )
2019-03-06 20:27:42 +05:30
real ( pReal ) , intent ( in ) :: a
real ( pReal ) , intent ( in ) , optional :: tol
real ( pReal ) :: eps
2019-03-07 23:59:48 +05:30
2019-03-06 20:27:42 +05:30
if ( present ( tol ) ) then
eps = tol
else
2019-03-07 15:32:27 +05:30
eps = PREAL_MIN * 1 0.0_pReal
2019-03-06 20:27:42 +05:30
endif
2019-09-23 21:59:38 +05:30
dEq0 = merge ( . True . , . False . , abs ( a ) < = eps )
2016-06-30 14:00:40 +05:30
end function dEq0
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Test floating point number with double precision for inequality to 0.
2017-11-21 13:48:03 +05:30
! replaces "!=0" but everything not representable as a normal number is treated as 0. Counterpart to dEq0
! https://de.mathworks.com/help/matlab/ref/realmin.html
! https://docs.oracle.com/cd/E19957-01/806-3568/ncg_math.html
2016-06-30 14:00:40 +05:30
!--------------------------------------------------------------------------------------------------
logical elemental pure function dNeq0 ( a , tol )
2019-03-06 20:27:42 +05:30
real ( pReal ) , intent ( in ) :: a
real ( pReal ) , intent ( in ) , optional :: tol
2019-03-07 23:59:48 +05:30
2019-03-06 20:27:42 +05:30
if ( present ( tol ) ) then
2019-09-23 12:23:56 +05:30
dNeq0 = . not . dEq0 ( a , tol )
2019-03-06 20:27:42 +05:30
else
2019-09-23 12:23:56 +05:30
dNeq0 = . not . dEq0 ( a )
2019-03-06 20:27:42 +05:30
endif
2016-06-30 14:00:40 +05:30
end function dNeq0
2016-05-29 14:15:03 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Test complex floating point numbers with double precision for equality.
2016-05-29 14:15:03 +05:30
! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq
2016-10-18 23:20:05 +05:30
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
2016-05-29 14:15:03 +05:30
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cEq ( a , b , tol )
2019-03-06 20:27:42 +05:30
complex ( pReal ) , intent ( in ) :: a , b
real ( pReal ) , intent ( in ) , optional :: tol
real ( pReal ) :: eps
2019-03-07 23:59:48 +05:30
2019-03-06 20:27:42 +05:30
if ( present ( tol ) ) then
eps = tol
else
2019-03-07 15:32:27 +05:30
eps = PREAL_EPSILON * maxval ( abs ( [ a , b ] ) )
2019-03-06 20:27:42 +05:30
endif
2019-09-23 21:59:38 +05:30
cEq = merge ( . True . , . False . , abs ( a - b ) < = eps )
2016-05-29 14:15:03 +05:30
end function cEq
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Test complex floating point numbers with double precision for inequality.
2016-05-29 14:15:03 +05:30
! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq
2016-10-18 23:20:05 +05:30
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
2016-05-29 14:15:03 +05:30
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cNeq ( a , b , tol )
2019-03-06 20:27:42 +05:30
complex ( pReal ) , intent ( in ) :: a , b
real ( pReal ) , intent ( in ) , optional :: tol
2019-03-07 23:59:48 +05:30
2019-03-06 20:27:42 +05:30
if ( present ( tol ) ) then
2019-09-23 12:23:56 +05:30
cNeq = . not . cEq ( a , b , tol )
2019-03-06 20:27:42 +05:30
else
2019-09-23 12:23:56 +05:30
cNeq = . not . cEq ( a , b )
2019-03-06 20:27:42 +05:30
endif
2016-05-29 14:15:03 +05:30
end function cNeq
2019-09-23 12:23:56 +05:30
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Decode byte array (C_SIGNED_CHAR) as C_FLOAT array (4 byte float).
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-13 23:44:34 +05:30
pure function prec_bytesToC_FLOAT ( bytes )
2020-09-06 21:06:05 +05:30
integer ( C_SIGNED_CHAR ) , dimension ( : ) , intent ( in ) :: bytes !< byte-wise representation of a C_FLOAT array
2020-09-12 12:44:33 +05:30
real ( C_FLOAT ) , dimension ( size ( bytes , kind = pI64 ) / ( storage_size ( 0._C_FLOAT , pI64 ) / 8_pI64 ) ) :: &
2020-09-13 23:44:34 +05:30
prec_bytesToC_FLOAT
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
prec_bytesToC_FLOAT = transfer ( bytes , prec_bytesToC_FLOAT , size ( prec_bytesToC_FLOAT ) )
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
end function prec_bytesToC_FLOAT
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Decode byte array (C_SIGNED_CHAR) as C_DOUBLE array (8 byte float).
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-13 23:44:34 +05:30
pure function prec_bytesToC_DOUBLE ( bytes )
2020-09-06 21:06:05 +05:30
integer ( C_SIGNED_CHAR ) , dimension ( : ) , intent ( in ) :: bytes !< byte-wise representation of a C_DOUBLE array
2020-09-12 12:44:33 +05:30
real ( C_DOUBLE ) , dimension ( size ( bytes , kind = pI64 ) / ( storage_size ( 0._C_DOUBLE , pI64 ) / 8_pI64 ) ) :: &
2020-09-13 23:44:34 +05:30
prec_bytesToC_DOUBLE
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
prec_bytesToC_DOUBLE = transfer ( bytes , prec_bytesToC_DOUBLE , size ( prec_bytesToC_DOUBLE ) )
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
end function prec_bytesToC_DOUBLE
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Decode byte array (C_SIGNED_CHAR) as C_INT32_T array (4 byte signed integer).
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-13 23:44:34 +05:30
pure function prec_bytesToC_INT32_T ( bytes )
2020-09-06 21:06:05 +05:30
integer ( C_SIGNED_CHAR ) , dimension ( : ) , intent ( in ) :: bytes !< byte-wise representation of a C_INT32_T array
2020-09-12 12:44:33 +05:30
integer ( C_INT32_T ) , dimension ( size ( bytes , kind = pI64 ) / ( storage_size ( 0_C_INT32_T , pI64 ) / 8_pI64 ) ) :: &
2020-09-13 23:44:34 +05:30
prec_bytesToC_INT32_T
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
prec_bytesToC_INT32_T = transfer ( bytes , prec_bytesToC_INT32_T , size ( prec_bytesToC_INT32_T ) )
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
end function prec_bytesToC_INT32_T
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Decode byte array (C_SIGNED_CHAR) as C_INT64_T array (8 byte signed integer).
2020-09-06 21:06:05 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-13 23:44:34 +05:30
pure function prec_bytesToC_INT64_T ( bytes )
2020-09-06 21:06:05 +05:30
integer ( C_SIGNED_CHAR ) , dimension ( : ) , intent ( in ) :: bytes !< byte-wise representation of a C_INT64_T array
2020-09-12 12:44:33 +05:30
integer ( C_INT64_T ) , dimension ( size ( bytes , kind = pI64 ) / ( storage_size ( 0_C_INT64_T , pI64 ) / 8_pI64 ) ) :: &
2020-09-13 23:44:34 +05:30
prec_bytesToC_INT64_T
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
prec_bytesToC_INT64_T = transfer ( bytes , prec_bytesToC_INT64_T , size ( prec_bytesToC_INT64_T ) )
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
end function prec_bytesToC_INT64_T
2020-09-06 21:06:05 +05:30
2019-09-23 12:23:56 +05:30
!--------------------------------------------------------------------------------------------------
2020-09-08 10:25:04 +05:30
!> @brief Check correctness of some prec functions.
2019-09-23 12:23:56 +05:30
!--------------------------------------------------------------------------------------------------
2020-05-16 20:35:03 +05:30
subroutine selfTest
2019-09-23 12:23:56 +05:30
integer , allocatable , dimension ( : ) :: realloc_lhs_test
2020-09-06 21:06:05 +05:30
real ( pReal ) , dimension ( 1 ) :: f
integer ( pInt ) , dimension ( 1 ) :: i
real ( pReal ) , dimension ( 2 ) :: r
2020-03-27 02:16:28 +05:30
2020-09-06 21:06:05 +05:30
realloc_lhs_test = [ 1 , 2 ]
2020-09-13 15:39:32 +05:30
if ( any ( realloc_lhs_test / = [ 1 , 2 ] ) ) error stop 'LHS allocation'
2020-09-06 21:06:05 +05:30
2019-09-23 12:23:56 +05:30
call random_number ( r )
r = r / minval ( r )
2020-09-13 15:39:32 +05:30
if ( . not . all ( dEq ( r , r + PREAL_EPSILON ) ) ) error stop 'dEq'
if ( dEq ( r ( 1 ) , r ( 2 ) ) . and . dNeq ( r ( 1 ) , r ( 2 ) ) ) error stop 'dNeq'
if ( . not . all ( dEq0 ( r - ( r + PREAL_MIN ) ) ) ) error stop 'dEq0'
2020-03-27 02:16:28 +05:30
2020-09-06 21:06:05 +05:30
! https://www.binaryconvert.com
! https://www.rapidtables.com/convert/number/binary-to-decimal.html
2020-09-13 23:44:34 +05:30
f = real ( prec_bytesToC_FLOAT ( int ( [ - 65 , + 11 , - 102 , + 75 ] , C_SIGNED_CHAR ) ) , pReal )
if ( dNeq ( f ( 1 ) , 2019110 2.0_pReal , 0.0_pReal ) ) error stop 'prec_bytesToC_FLOAT'
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
f = real ( prec_bytesToC_DOUBLE ( int ( [ 0 , 0 , 0 , - 32 , + 119 , + 65 , + 115 , 65 ] , C_SIGNED_CHAR ) ) , pReal )
if ( dNeq ( f ( 1 ) , 2019110 2.0_pReal , 0.0_pReal ) ) error stop 'prec_bytesToC_DOUBLE'
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
i = int ( prec_bytesToC_INT32_T ( int ( [ + 126 , + 23 , + 52 , + 1 ] , C_SIGNED_CHAR ) ) , pInt )
if ( i ( 1 ) / = 20191102_pInt ) error stop 'prec_bytesToC_INT32_T'
2020-09-06 21:06:05 +05:30
2020-09-13 23:44:34 +05:30
i = int ( prec_bytesToC_INT64_T ( int ( [ + 126 , + 23 , + 52 , + 1 , 0 , 0 , 0 , 0 ] , C_SIGNED_CHAR ) ) , pInt )
if ( i ( 1 ) / = 20191102_pInt ) error stop 'prec_bytesToC_INT64_T'
2020-03-27 02:16:28 +05:30
2020-05-16 20:35:03 +05:30
end subroutine selfTest
2019-09-23 12:23:56 +05:30
2012-03-06 20:22:48 +05:30
end module prec