prec.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/prec.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/prec.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module prec
  use, intrinsic :: ieee_arithmetic
 
  implicit none
  public
 
  ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
  integer,     parameter :: preal      = ieee_selected_real_kind(15,307)                            
 
 
 
  integer,     parameter :: pint       = selected_int_kind(9)                                       
 
  integer,     parameter :: plongint   = selected_int_kind(18)                                      
  integer,     parameter :: pstringlen = 256
  integer,     parameter :: ppathlen   = 4096
 
  real(preal), parameter :: tol_math_check = 1.0e-8_preal                                           
 
 
  type :: group_float                                                                               
    real(preal), dimension(:), pointer :: p
  end type group_float
 
  type :: group_int
    integer, dimension(:), pointer :: p
  end type group_int
 
  ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
  type :: tstate
    integer :: &
      sizestate        = 0, &                                                                       !< size of state
      sizedotstate     = 0, &                                                                       
      offsetdeltastate = 0, &                                                                       
      sizedeltastate   = 0                                                                          
    real(preal), pointer,     dimension(:), contiguous :: &
      atol
    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
      state0, &
      state, &                                                                                      !< state
      dotstate, &                                                                                   !< rate of state change
      deltastate
    real(preal), allocatable, dimension(:,:) :: &
      partionedstate0, &
      substate0, &
      previousdotstate, &
      previousdotstate2
    real(preal), allocatable, dimension(:,:,:) :: &
      rk4dotstate, &
      rkck45dotstate
  end type
 
  type, extends(tstate) :: tplasticstate
    logical :: &
      nonlocal = .false.
    real(preal), pointer,     dimension(:,:) :: &
      sliprate, &                                                                                   !< slip rate
      accumulatedslip
  end type
 
  type :: tsourcestate
    type(tstate), dimension(:), allocatable :: p
  end type
 
  type :: thomogmapping
    integer, pointer, dimension(:,:) :: p
  end type
 
  real(preal), private, parameter :: preal_epsilon = epsilon(0.0_preal)                             
  real(preal), private, parameter :: preal_min     = tiny(0.0_preal)                                
 
  integer,                   dimension(0), parameter :: &
    emptyintarray    = [integer::]
  real(preal),               dimension(0), parameter :: &
    emptyrealarray   = [real(preal)::]
  character(len=pStringLen), dimension(0), parameter :: &
    emptystringarray = [character(len=pstringlen)::]
 
  private :: &
    unittest
 
contains
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine prec_init
 
  write(6,'(/,a)') ' <<<+-  prec init  -+>>>'
 
  write(6,'(a,i3)')    ' Size of integer in bit: ',bit_size(0)
  write(6,'(a,i19)')   '   Maximum value:        ',huge(0)
  write(6,'(/,a,i3)')  ' Size of float in bit:   ',storage_size(0.0_preal)
  write(6,'(a,e10.3)') '   Maximum value:        ',huge(0.0_preal)
  write(6,'(a,e10.3)') '   Minimum value:        ',tiny(0.0_preal)
  write(6,'(a,i3)')    '   Decimal precision:    ',precision(0.0_preal)
 
  call unittest
 
end subroutine prec_init
 
 
!--------------------------------------------------------------------------------------------------
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! AlmostEqualRelative
!--------------------------------------------------------------------------------------------------
logical elemental pure function dEq(a,b,tol)
 
  real(preal), intent(in)           :: a,b
  real(preal), intent(in), optional :: tol
  real(preal)                       :: eps
 
  if (present(tol)) then
    eps = tol
  else
    eps = preal_epsilon * maxval(abs([a,b]))
  endif
 
  deq = merge(.true.,.false.,abs(a-b) <= eps)
 
end function deq
 
 
!--------------------------------------------------------------------------------------------------
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! AlmostEqualRelative NOT
!--------------------------------------------------------------------------------------------------
logical elemental pure function dneq(a,b,tol)
 
  real(preal), intent(in)           :: a,b
  real(preal), intent(in), optional :: tol
 
  if (present(tol)) then
    dneq = .not. deq(a,b,tol)
  else
    dneq = .not. deq(a,b)
  endif
 
end function dneq
 
 
!--------------------------------------------------------------------------------------------------
! 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
!--------------------------------------------------------------------------------------------------
logical elemental pure function deq0(a,tol)
 
  real(preal), intent(in)           :: a
  real(preal), intent(in), optional :: tol
  real(preal)                       :: eps
 
  if (present(tol)) then
    eps = tol
  else
    eps = preal_min * 10.0_preal
  endif
 
  deq0 = merge(.true.,.false.,abs(a) <= eps)
 
end function deq0
 
 
!--------------------------------------------------------------------------------------------------
! 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
!--------------------------------------------------------------------------------------------------
logical elemental pure function dneq0(a,tol)
 
  real(preal), intent(in)           :: a
  real(preal), intent(in), optional :: tol
 
  if (present(tol)) then
    dneq0 = .not. deq0(a,tol)
  else
    dneq0 = .not. deq0(a)
  endif
 
end function dneq0
 
 
!--------------------------------------------------------------------------------------------------
! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function ceq(a,b,tol)
 
  complex(pReal), intent(in)           :: a,b
  real(preal),    intent(in), optional :: tol
  real(preal)                          :: eps
 
  if (present(tol)) then
    eps = tol
  else
    eps = preal_epsilon * maxval(abs([a,b]))
  endif
 
  ceq = merge(.true.,.false.,abs(a-b) <= eps)
 
end function ceq
 
 
!--------------------------------------------------------------------------------------------------
! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq
! https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
! probably a component wise comparison would be more accurate than the comparsion of the absolute
! value
!--------------------------------------------------------------------------------------------------
logical elemental pure function cneq(a,b,tol)
 
  complex(pReal), intent(in)           :: a,b
  real(preal),    intent(in), optional :: tol
 
  if (present(tol)) then
    cneq = .not. ceq(a,b,tol)
  else
    cneq = .not. ceq(a,b)
  endif
 
end function cneq
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine unittest
 
  integer, allocatable, dimension(:) :: realloc_lhs_test
  real(preal), dimension(2) :: r
  external :: &
    quit
 
  call random_number(r)
  r = r/minval(r)
  if(.not. all(deq(r,r+preal_epsilon)))    call quit(9000)
  if(deq(r(1),r(2)) .and. dneq(r(1),r(2))) call quit(9000)
  if(.not. all(deq0(r-(r+preal_min))))     call quit(9000)
 
  realloc_lhs_test = [1,2]
  if (any(realloc_lhs_test/=[1,2])) call quit(9000)
 
end subroutine unittest
 
end module prec