223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

added error for to IO to substitute stop statement in kdtree2 

explicitly defined all functions in as either public or private in the modules to have a quick overview on all functions and parameters that are available
This commit is contained in:
Martin Diehl 2012-03-06 14:52:48 +00:00
parent 9b17015b5a
commit d00c3c9e19
6 changed files with 2464 additions and 2366 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

188
code/DAMASK_spectral_interface.f90
View File

@ -1,7 +1,7 @@
! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH ! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
! !
! This file is part of DAMASK, ! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit. ! the Düsseldorf Advanced Material Simulation Kit.
! !
! DAMASK is free software: you can redistribute it and/or modify ! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by ! it under the terms of the GNU General Public License as published by
@ -16,39 +16,60 @@
! You should have received a copy of the GNU General Public License ! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>. ! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
! !
!############################################################## !--------------------------------------------------------------------------------------------------
!* $Id$ !* $Id$
!******************************************************************** !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Interfacing between the spectral solver and the material subroutines provided
!! by DAMASK
!--------------------------------------------------------------------------------------------------
module DAMASK_interface
MODULE DAMASK_interface
implicit none implicit none
private
character(len=64), parameter, public :: FEsolver = 'Spectral' !> Keyword for spectral solver
character(len=5), parameter, public :: inputFileExtension = '.geom' !> File extension for geometry description
character(len=4), parameter, public :: logFileExtension = '.log' !> Dummy variable as the spectral solver has no log
character(len=1024), private :: geometryParameter, &
loadcaseParameter
character(len=64), parameter :: FEsolver = 'Spectral' public :: getSolverWorkingDirectoryName, &
character(len=5), parameter :: InputFileExtension = '.geom' getSolverJobName, &
character(len=4), parameter :: LogFileExtension = '.log' !until now, we don't have a log file. But IO.f90 requires it getLoadCase, &
character(len=1024) :: geometryParameter,loadcaseParameter getLoadCaseName, &
CONTAINS getModelName, &
DAMASK_interface_init
private :: rectifyPath, &
makeRelativePath, &
getPathSep
!******************************************************************** contains
! initialize interface module
! !--------------------------------------------------------------------------------------------------
!******************************************************************** !> @brief Initializes the solver by interpreting the command line arguments. Also writes
subroutine DAMASK_interface_init() !! information on computation on screen
!--------------------------------------------------------------------------------------------------
subroutine DAMASK_interface_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt use prec, only: pInt
implicit none
character(len=1024) commandLine, hostName, userName implicit none
integer :: i, start = 0, length=0 character(len=1024) :: commandLine, & !> command line call as string
integer, dimension(8) :: date_and_time_values ! type default integer hostName, & !> name of computer
userName !> name of user calling the executable
integer :: i, &
start = 0,&
length=0
integer, dimension(8) :: dateAndTime ! type default integer
call get_command(commandLine) call get_command(commandLine)
call DATE_AND_TIME(VALUES=date_and_time_values) call date_and_time(values = dateAndTime)
do i=1,len(commandLine) ! remove capitals do i = 1,len(commandLine) ! remove capitals
if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) commandLine(i:i) = & if(64<iachar(commandLine(i:i)) .and. iachar(commandLine(i:i))<91) &
achar(iachar(commandLine(i:i))+32) commandLine(i:i) = achar(iachar(commandLine(i:i))+32)
enddo enddo
if(index(commandLine,' -h ',.true.)>0 .or. index(commandLine,' --help ',.true.)>0) then ! search for ' -h ' or '--help' if(index(commandLine,' -h ',.true.) > 0 .or. index(commandLine,' --help ',.true.) > 0) then ! search for ' -h ' or '--help'
write(6,*) '$Id$' write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
print '(a)', '#############################################################' print '(a)', '#############################################################'
@ -146,12 +167,12 @@ subroutine DAMASK_interface_init()
write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>' write(6,*) '<<<+- DAMASK_spectral_interface init -+>>>'
write(6,*) '$Id$' write(6,*) '$Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',date_and_time_values(3),'/',& write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',&
date_and_time_values(2),'/',& dateAndTime(2),'/',&
date_and_time_values(1) dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',date_and_time_values(5),':',& write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',&
date_and_time_values(6),':',& dateAndTime(6),':',&
date_and_time_values(7) dateAndTime(7)
write(6,*) 'Host Name: ', trim(hostName) write(6,*) 'Host Name: ', trim(hostName)
write(6,*) 'User Name: ', trim(userName) write(6,*) 'User Name: ', trim(userName)
write(6,*) 'Path Separator: ', getPathSep() write(6,*) 'Path Separator: ', getPathSep()
@ -160,17 +181,15 @@ subroutine DAMASK_interface_init()
write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter) write(6,*) 'Loadcase Parameter: ', trim(loadcaseParameter)
if (start/=3_pInt) write(6,*) 'Restart Parameter: ', trim(commandLine(start:start+length)) if (start/=3_pInt) write(6,*) 'Restart Parameter: ', trim(commandLine(start:start+length))
endsubroutine DAMASK_interface_init end subroutine DAMASK_interface_init
!******************************************************************** !--------------------------------------------------------------------------------------------------
! extract working directory from loadcase file !> @brief extract working directory from loadcase file possibly based on current working dir
! possibly based on current working dir !--------------------------------------------------------------------------------------------------
!********************************************************************
function getSolverWorkingDirectoryName() function getSolverWorkingDirectoryName()
implicit none implicit none
character(len=1024) :: cwd, getSolverWorkingDirectoryName
character(len=1024) cwd,getSolverWorkingDirectoryName
character :: pathSep character :: pathSep
pathSep = getPathSep() pathSep = getPathSep()
@ -184,33 +203,29 @@ function getSolverWorkingDirectoryName()
getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName) getSolverWorkingDirectoryName = rectifyPath(getSolverWorkingDirectoryName)
endfunction getSolverWorkingDirectoryName end function getSolverWorkingDirectoryName
!********************************************************************
! basename of geometry file from command line arguments !--------------------------------------------------------------------------------------------------
! !> @brief basename of geometry file from command line arguments
!******************************************************************** !--------------------------------------------------------------------------------------------------
function getSolverJobName() character(len=1024) function getSolverJobName()
implicit none implicit none
character(1024) :: getSolverJobName
getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase()) getSolverJobName = trim(getModelName())//'_'//trim(getLoadCase())
endfunction getSolverJobName end function getSolverJobName
!********************************************************************
! basename of geometry file from command line arguments !--------------------------------------------------------------------------------------------------
! !> @brief basename of geometry file from command line arguments
!******************************************************************** !--------------------------------------------------------------------------------------------------
function getModelName() character(len=1024) function getModelName()
use prec, only: pInt use prec, only: pInt
implicit none implicit none
character(len=1024) :: cwd
character(1024) getModelName, cwd
integer :: posExt,posSep integer :: posExt,posSep
character :: pathSep character :: pathSep
@ -231,17 +246,15 @@ function getModelName()
getModelName = makeRelativePath(getSolverWorkingDirectoryName(),& getModelName = makeRelativePath(getSolverWorkingDirectoryName(),&
getModelName) getModelName)
endfunction getModelName end function getModelName
!********************************************************************
! name of load case file exluding extension !--------------------------------------------------------------------------------------------------
! !> @brief name of load case file exluding extension
!******************************************************************** !--------------------------------------------------------------------------------------------------
function getLoadCase() character(len=1024) function getLoadCase()
implicit none implicit none
character(1024) :: getLoadCase
integer :: posExt,posSep integer :: posExt,posSep
character :: pathSep character :: pathSep
@ -252,18 +265,16 @@ function getLoadCase()
if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present if (posExt <= posSep) posExt = len_trim(loadcaseParameter)+1 ! no extension present
getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension getLoadCase = loadcaseParameter(posSep+1:posExt-1) ! name of load case file exluding extension
endfunction getLoadCase end function getLoadCase
!******************************************************************** !--------------------------------------------------------------------------------------------------
! relative path of loadcase from command line arguments !> @brief relative path of loadcase from command line arguments
! !--------------------------------------------------------------------------------------------------
!******************************************************************** character(len=1024) function getLoadcaseName()
function getLoadcaseName()
implicit none implicit none
character(len=1024) :: cwd
character(len=1024) :: getLoadcaseName,cwd
integer :: posExt = 0, posSep integer :: posExt = 0, posSep
character :: pathSep character :: pathSep
@ -282,17 +293,15 @@ function getLoadcaseName()
getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),& getLoadcaseName = makeRelativePath(getSolverWorkingDirectoryName(),&
getLoadcaseName) getLoadcaseName)
endfunction getLoadcaseName end function getLoadcaseName
!******************************************************************** !--------------------------------------------------------------------------------------------------
! remove ../ and ./ from path !> @brief remove ../ and ./ from path
! !--------------------------------------------------------------------------------------------------
!********************************************************************
function rectifyPath(path) function rectifyPath(path)
implicit none implicit none
character(len=*) :: path character(len=*) :: path
character(len=len_trim(path)) :: rectifyPath character(len=len_trim(path)) :: rectifyPath
character :: pathSep character :: pathSep
@ -324,19 +333,16 @@ function rectifyPath(path)
enddo enddo
if(len_trim(rectifyPath) == 0) rectifyPath = pathSep if(len_trim(rectifyPath) == 0) rectifyPath = pathSep
end function rectifyPath end function rectifyPath
!******************************************************************** !--------------------------------------------------------------------------------------------------
! relative path from absolute a to absolute b !> @brief relative path from absolute a to absolute b
! !--------------------------------------------------------------------------------------------------
!******************************************************************** character(len=1024) function makeRelativePath(a,b)
function makeRelativePath(a,b)
implicit none implicit none
character (len=*) :: a,b character (len=*) :: a,b
character (len=1024) :: makeRelativePath
character :: pathSep character :: pathSep
integer :: i,posLastCommonSlash,remainingSlashes !no pInt integer :: i,posLastCommonSlash,remainingSlashes !no pInt
@ -353,18 +359,18 @@ function makeRelativePath(a,b)
enddo enddo
makeRelativePath = repeat('..'//pathSep,remainingSlashes)//b(posLastCommonSlash+1:len_trim(b)) makeRelativePath = repeat('..'//pathSep,remainingSlashes)//b(posLastCommonSlash+1:len_trim(b))
endfunction makeRelativePath end function makeRelativePath
!******************************************************************** !--------------------------------------------------------------------------------------------------
! counting / and \ in $PATH System variable !> @brief counting / and \ in $PATH System variable the character occuring more often is assumed
! the character occuring more often is assumed to be the path separator !! to be the path separator
!******************************************************************** !--------------------------------------------------------------------------------------------------
function getPathSep() character function getPathSep()
use prec, only: pInt use prec, only: pInt
implicit none implicit none
character :: getPathSep
character(len=2048) path character(len=2048) path
integer(pInt) :: backslash = 0_pInt, slash = 0_pInt integer(pInt) :: backslash = 0_pInt, slash = 0_pInt
integer :: i integer :: i
@ -383,4 +389,4 @@ function getPathSep()
end function end function
END MODULE end module

1041
code/IO.f90
View File

File diff suppressed because it is too large Load Diff

25
code/Makefile
View File

@ -135,8 +135,7 @@ COMPILE_OPTIONS_ifort :=-fpp\
-warn ignore_loc\ -warn ignore_loc\
-warn alignments\ -warn alignments\
-warn unused\ -warn unused\
-warn errors\ -warn errors
-warn stderrors
endif endif
#-fpp: preprocessor #-fpp: preprocessor
@ -152,7 +151,7 @@ endif
# alignments: data that is not naturally aligned # alignments: data that is not naturally aligned
# unused: declared variables that are never used # unused: declared variables that are never used
# errors: warnings are changed to errors # errors: warnings are changed to errors
# stderrors: warnings about Fortran standard violations are changed to errors # stderrors: warnings about Fortran standard violations are changed to errors (STANDARD_CHECK)
# #
################################################################################################### ###################################################################################################
#MORE OPTIONS FOR DEBUGGING DURING COMPILING #MORE OPTIONS FOR DEBUGGING DURING COMPILING
@ -182,7 +181,8 @@ endif
################################################################################################### ###################################################################################################
ifeq "$(FASTBUILD)" "YES" ifeq "$(FASTBUILD)" "YES"
COMPILE_OPTIONS_gfortran :=-xf95-cpp-input COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
-fno-range-check
else else
COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\ COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
-ffree-line-length-132\ -ffree-line-length-132\
@ -205,7 +205,8 @@ COMPILE_OPTIONS_gfortran :=-xf95-cpp-input\
-Wunsafe-loop-optimizations\ -Wunsafe-loop-optimizations\
-Wunused\ -Wunused\
-Wall\ -Wall\
-Wextra -Wextra\
-Wintrinsics-std
endif endif
#-xf95-cpp-input: preprocessor #-xf95-cpp-input: preprocessor
@ -231,20 +232,20 @@ endif
# -value: # -value:
# -parameter: find usused variables with "parameter" attribute # -parameter: find usused variables with "parameter" attribute
#-Wextra: #-Wextra:
#-Wintrinsics-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers:
################################################################################################### ###################################################################################################
#OPTIONS FOR GFORTRAN 4.6 #OPTIONS FOR GFORTRAN 4.6
#-Wsuggest-attribute=const: #-Wsuggest-attribute=const:
#-Wsuggest-attribute=noreturn: #-Wsuggest-attribute=noreturn:
#-Wsuggest-attribute=pure: #-Wsuggest-attribute=pure:
#-Wreal-q-constant: Warn about real-literal-constants with 'q' exponent-letter #-Wreal-q-constant: Warn about real-literal-constants with 'q' exponent-letter
#
#MORE OPTIONS FOR DEBUGGING DURING COMPILING #MORE OPTIONS FOR DEBUGGING DURING COMPILING
#-Wline-truncation: too many warnings because we have comments beyond character 132 #-Wline-truncation: too many warnings because we have comments beyond character 132
#-Wintrinsic-std: warnings because of "flush" is not longer in the standard, but still an intrinsic fuction of the compilers: #-Warray-temporarieswarnings: because we have many temporary arrays (performance issue?):
#-Warray-temporarieswarnings: #-Wimplicit-interface:
# because we have many temporary arrays (performance issue?): #-pedantic-errors:
#-Wimplicit-interface #-fmodule-private:
#-pedantic-errors
#-fmodule-private
# #
#OPTIONS FOR DEGUBBING DURING RUNTIME #OPTIONS FOR DEGUBBING DURING RUNTIME
#-fcheck-bounds: check if an array index is too small (<1) or too large! #-fcheck-bounds: check if an array index is too small (<1) or too large!
@ -264,7 +265,7 @@ COMPILED_FILES = prec.o DAMASK_spectral_interface.o IO.o numerics.o debug.o math
homogenization_RGC.o homogenization_isostrain.o homogenization.o CPFEM.o crystallite.o homogenization_RGC.o homogenization_isostrain.o homogenization.o CPFEM.o crystallite.o
DAMASK_spectral.exe: DAMASK_spectral.o DAMASK_spectral.exe: DAMASK_spectral.o
$(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90))$(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \ $(PREFIX) $(COMPILERNAME) $(OPENMP_FLAG_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(STANDARD_CHECK_$(F90)) \
-o DAMASK_spectral.exe DAMASK_spectral.o \ -o DAMASK_spectral.exe DAMASK_spectral.o \
$(COMPILED_FILES) $(LIB_DIRS) $(LIBRARIES) $(SUFFIX) $(COMPILED_FILES) $(LIB_DIRS) $(LIBRARIES) $(SUFFIX)

231
code/kdtree2.f90
View File

@ -12,12 +12,15 @@
!####################################################### !#######################################################
module kdtree2_priority_queue_module module kdtree2_priority_queue_module
use prec use prec, only: pInt, &
pReal
! !
! maintain a priority queue (PQ) of data, pairs of 'priority/payload', ! maintain a priority queue (PQ) of data, pairs of 'priority/payload',
! implemented with a binary heap. This is the type, and the 'dis' field ! implemented with a binary heap. This is the type, and the 'dis' field
! is the priority. ! is the priority.
! !
implicit none
private
type kdtree2_result type kdtree2_result
! a pair of distances, indexes ! a pair of distances, indexes
real(pReal) :: dis = 0.0_pReal real(pReal) :: dis = 0.0_pReal
@ -74,12 +77,10 @@ module kdtree2_priority_queue_module
public :: pq_create public :: pq_create
public :: pq_delete, pq_insert public :: pq_delete, pq_insert
public :: pq_extract_max, pq_max, pq_replace_max, pq_maxpri public :: pq_extract_max, pq_max, pq_replace_max, pq_maxpri
private
contains contains
function pq_create(results_in) result(res)
function pq_create(results_in) result(res)
! !
! Create a priority queue from ALREADY allocated ! Create a priority queue from ALREADY allocated
! array pointers for storage. NOTE! It will NOT ! array pointers for storage. NOTE! It will NOT
@ -93,6 +94,7 @@ contains
! allocate(x(1000),k(1000)) ! allocate(x(1000),k(1000))
! pq => pq_create(x,k) ! pq => pq_create(x,k)
! !
implicit none
type(kdtree2_result), target:: results_in(:) type(kdtree2_result), target:: results_in(:)
type(pq) :: res type(pq) :: res
! !
@ -106,7 +108,7 @@ contains
res%elems => results_in res%elems => results_in
res%heap_size = 0_pInt res%heap_size = 0_pInt
return return
end function pq_create end function pq_create
! !
! operations for getting parents and left + right children ! operations for getting parents and left + right children
@ -141,12 +143,13 @@ contains
! return ! return
! end function compare_priority ! end function compare_priority
subroutine heapify(a,i_in) subroutine heapify(a,i_in)
! !
! take a heap rooted at 'i' and force it to be in the ! take a heap rooted at 'i' and force it to be in the
! heap canonical form. This is performance critical ! heap canonical form. This is performance critical
! and has been tweaked a little to reflect this. ! and has been tweaked a little to reflect this.
! !
implicit none
type(pq),pointer :: a type(pq),pointer :: a
integer(pInt), intent(in) :: i_in integer(pInt), intent(in) :: i_in
! !
@ -159,7 +162,7 @@ contains
i = i_in i = i_in
bigloop: do bigloop: do
l = 2_pInt*i ! left(i) l = 2_pInt*i ! left(i)
r = l+1_pInt ! right(i) r = l+1_pInt ! right(i)
! !
@ -216,44 +219,51 @@ bigloop: do
end if end if
enddo bigloop enddo bigloop
return return
end subroutine heapify end subroutine heapify
subroutine pq_max(a,e) subroutine pq_max(a,e)
! !
! return the priority and its payload of the maximum priority element ! return the priority and its payload of the maximum priority element
! on the queue, which should be the first one, if it is ! on the queue, which should be the first one, if it is
! in heapified form. ! in heapified form.
! !
use IO, only: IO_error
implicit none
type(pq),pointer :: a type(pq),pointer :: a
type(kdtree2_result),intent(out) :: e type(kdtree2_result),intent(out) :: e
if (a%heap_size .gt. 0) then if (a%heap_size .gt. 0) then
e = a%elems(1) e = a%elems(1)
else else
write (*,*) 'PQ_MAX: ERROR, heap_size < 1' call IO_error (500_pInt, ext_msg='PQ_MAX: heap_size < 1')
stop
endif endif
return return
end subroutine pq_max end subroutine pq_max
real(pReal) function pq_maxpri(a) real(pReal) function pq_maxpri(a)
use IO, only: IO_error
implicit none
type(pq), pointer :: a type(pq), pointer :: a
if (a%heap_size .gt. 0) then if (a%heap_size .gt. 0) then
pq_maxpri = a%elems(1)%dis pq_maxpri = a%elems(1)%dis
else else
write (*,*) 'PQ_MAX_PRI: ERROR, heapsize < 1' call IO_error (500_pInt,ext_msg='PPQ_MAX_PRI: heap_size < 1')
stop
endif endif
return return
end function pq_maxpri end function pq_maxpri
subroutine pq_extract_max(a,e) subroutine pq_extract_max(a,e)
! !
! return the priority and payload of maximum priority ! return the priority and payload of maximum priority
! element, and remove it from the queue. ! element, and remove it from the queue.
! (equivalent to 'pop()' on a stack) ! (equivalent to 'pop()' on a stack)
! !
use IO, only: IO_error
implicit none
type(pq),pointer :: a type(pq),pointer :: a
type(kdtree2_result), intent(out) :: e type(kdtree2_result), intent(out) :: e
@ -271,18 +281,17 @@ bigloop: do
call heapify(a,1_pInt) call heapify(a,1_pInt)
return return
else else
write (*,*) 'PQ_EXTRACT_MAX: error, attempted to pop non-positive PQ' call IO_error (500_pInt,ext_msg='PQ_EXTRACT_MAX: attempted to pop non-positive PQ')
stop
end if end if
end subroutine pq_extract_max
end subroutine pq_extract_max
real(pReal) function pq_insert(a,dis,idx) real(pReal) function pq_insert(a,dis,idx)
! !
! Insert a new element and return the new maximum priority, ! Insert a new element and return the new maximum priority,
! which may or may not be the same as the old maximum priority. ! which may or may not be the same as the old maximum priority.
! !
implicit none
type(pq),pointer :: a type(pq),pointer :: a
real(pReal), intent(in) :: dis real(pReal), intent(in) :: dis
integer(pInt), intent(in) :: idx integer(pInt), intent(in) :: idx
@ -320,10 +329,11 @@ bigloop: do
return return
! end if ! end if
end function pq_insert end function pq_insert
subroutine pq_adjust_heap(a,i) subroutine pq_adjust_heap(a,i)
type(pq),pointer :: a implicit none
type(pq), pointer :: a
integer(pInt), intent(in) :: i integer(pInt), intent(in) :: i
! !
! nominally arguments (a,i), but specialize for a=1 ! nominally arguments (a,i), but specialize for a=1
@ -360,20 +370,21 @@ bigloop: do
end do end do
a%elems(parent) = e a%elems(parent) = e
return return
end subroutine pq_adjust_heap end subroutine pq_adjust_heap
real(pReal) function pq_replace_max(a,dis,idx) real(pReal) function pq_replace_max(a,dis,idx)
! !
! Replace the extant maximum priority element ! Replace the extant maximum priority element
! in the PQ with (dis,idx). Return ! in the PQ with (dis,idx). Return
! the new maximum priority, which may be larger ! the new maximum priority, which may be larger
! or smaller than the old one. ! or smaller than the old one.
! !
implicit none
type(pq),pointer :: a type(pq),pointer :: a
real(pReal), intent(in) :: dis real(pReal), intent(in) :: dis
integer(pInt), intent(in) :: idx integer(pInt), intent(in) :: idx
! type(kdtree2_result), intent(in) :: e !type(kdtree2_result), intent(in) :: e
! not tested as well! ! not tested as well!
integer(pInt) :: parent, child, N integer(pInt) :: parent, child, N
@ -432,18 +443,20 @@ bigloop: do
pq_replace_max = pq_insert(a,dis,idx) pq_replace_max = pq_insert(a,dis,idx)
endif endif
return return
end function pq_replace_max end function pq_replace_max
subroutine pq_delete(a,i) subroutine pq_delete(a,i)
! !
! delete item with index 'i' ! delete item with index 'i'
! !
use IO, only: IO_error
implicit none
type(pq),pointer :: a type(pq),pointer :: a
integer(pInt) :: i integer(pInt) :: i
if ((i .lt. 1) .or. (i .gt. a%heap_size)) then if ((i .lt. 1) .or. (i .gt. a%heap_size)) then
write (*,*) 'PQ_DELETE: error, attempt to remove out of bounds element.' call IO_error (500_pInt,ext_msg='PQ_DELETE: attempt to remove out of bounds element')
stop
endif endif
! swap the item to be deleted with the last element ! swap the item to be deleted with the last element
@ -453,7 +466,7 @@ bigloop: do
call heapify(a,i) call heapify(a,i)
end subroutine pq_delete end subroutine pq_delete
end module kdtree2_priority_queue_module end module kdtree2_priority_queue_module
@ -461,6 +474,9 @@ end module kdtree2_priority_queue_module
module kdtree2_module module kdtree2_module
use prec use prec
use kdtree2_priority_queue_module use kdtree2_priority_queue_module
implicit none
private
! K-D tree routines in Fortran 90 by Matt Kennel. ! K-D tree routines in Fortran 90 by Matt Kennel.
! Original program was written in Sather by Steve Omohundro and ! Original program was written in Sather by Steve Omohundro and
! Matt Kennel. Only the Euclidean metric is supported. ! Matt Kennel. Only the Euclidean metric is supported.
@ -583,14 +599,11 @@ module kdtree2_module
integer(pInt), pointer :: ind(:) ! temp pointer to indexes integer(pInt), pointer :: ind(:) ! temp pointer to indexes
end type tree_search_record end type tree_search_record
private
! everything else is private.
type(tree_search_record), save, target :: sr ! A GLOBAL VARIABLE for search type(tree_search_record), save, target :: sr ! A GLOBAL VARIABLE for search
contains contains
function kdtree2_create(input_data,myDim,sort,rearrange) result (mr) function kdtree2_create(input_data,myDim,sort,rearrange) result (mr)
! !
! create the actual tree structure, given an input array of data. ! create the actual tree structure, given an input array of data.
! !
@ -614,6 +627,9 @@ contains
! building takes longer, and extra memory is used. ! building takes longer, and extra memory is used.
! !
! .. Function Return Cut_value .. ! .. Function Return Cut_value ..
use IO, only: IO_error
implicit none
type (kdtree2), pointer :: mr type (kdtree2), pointer :: mr
integer(pInt), intent(in), optional :: myDim integer(pInt), intent(in), optional :: myDim
logical, intent(in), optional :: sort logical, intent(in), optional :: sort
@ -643,7 +659,7 @@ contains
write (*,*) 'KD_TREE_TRANS: note, that new format is myData(1:D,1:N)' write (*,*) 'KD_TREE_TRANS: note, that new format is myData(1:D,1:N)'
write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree' write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree'
write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.'
stop call IO_error (500_pInt)
end if end if
call build_tree(mr) call build_tree(mr)
@ -670,9 +686,10 @@ contains
nullify(mr%rearranged_data) nullify(mr%rearranged_data)
endif endif
end function kdtree2_create end function kdtree2_create
subroutine build_tree(tp) subroutine build_tree(tp)
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
! .. ! ..
integer(pInt) :: j integer(pInt) :: j
@ -683,10 +700,11 @@ contains
tp%ind(j) = j tp%ind(j) = j
end forall end forall
tp%root => build_tree_for_range(tp,1_pInt,tp%n, dummy) tp%root => build_tree_for_range(tp,1_pInt,tp%n, dummy)
end subroutine build_tree end subroutine build_tree
recursive function build_tree_for_range(tp,l,u,parent) result (res) recursive function build_tree_for_range(tp,l,u,parent) result (res)
! .. Function Return Cut_value .. ! .. Function Return Cut_value ..
implicit none
type (tree_node), pointer :: res type (tree_node), pointer :: res
! .. ! ..
! .. Structure Arguments .. ! .. Structure Arguments ..
@ -821,10 +839,9 @@ contains
res%box%lower = min(res%left%box%lower,res%right%box%lower) res%box%lower = min(res%left%box%lower,res%right%box%lower)
endif endif
end if end if
end function build_tree_for_range end function build_tree_for_range
integer(pInt) function select_on_coordinate_value(v,ind,c,alpha,li,ui) & integer(pInt) function select_on_coordinate_value(v,ind,c,alpha,li,ui) result(res)
result(res)
! Move elts of ind around between l and u, so that all points ! Move elts of ind around between l and u, so that all points
! <= than alpha (in c cooordinate) are first, and then ! <= than alpha (in c cooordinate) are first, and then
! all points > alpha are second. ! all points > alpha are second.
@ -842,6 +859,7 @@ contains
! The algorithm finishes when the unknown stack is empty. ! The algorithm finishes when the unknown stack is empty.
! !
! .. Scalar Arguments .. ! .. Scalar Arguments ..
implicit none
integer(pInt), intent (In) :: c, li, ui integer(pInt), intent (In) :: c, li, ui
real(pReal), intent(in) :: alpha real(pReal), intent(in) :: alpha
! .. ! ..
@ -882,13 +900,14 @@ contains
res = lb-1_pInt res = lb-1_pInt
endif endif
end function select_on_coordinate_value end function select_on_coordinate_value
subroutine select_on_coordinate(v,ind,c,k,li,ui) subroutine select_on_coordinate(v,ind,c,k,li,ui)
! Move elts of ind around between l and u, so that the kth ! Move elts of ind around between l and u, so that the kth
! element ! element
! is >= those below, <= those above, in the coordinate c. ! is >= those below, <= those above, in the coordinate c.
! .. Scalar Arguments .. ! .. Scalar Arguments ..
implicit none
integer(pInt), intent (In) :: c, k, li, ui integer(pInt), intent (In) :: c, k, li, ui
! .. ! ..
integer(pInt) :: i, l, m, s, t, u integer(pInt) :: i, l, m, s, t, u
@ -915,14 +934,15 @@ contains
if (m<=k) l = m + 1_pInt if (m<=k) l = m + 1_pInt
if (m>=k) u = m - 1_pInt if (m>=k) u = m - 1_pInt
end do end do
end subroutine select_on_coordinate end subroutine select_on_coordinate
subroutine spread_in_coordinate(tp,c,l,u,interv) subroutine spread_in_coordinate(tp,c,l,u,interv)
! the spread in coordinate 'c', between l and u. ! the spread in coordinate 'c', between l and u.
! !
! Return lower bound in 'smin', and upper in 'smax', ! Return lower bound in 'smin', and upper in 'smax',
! .. ! ..
! .. Structure Arguments .. ! .. Structure Arguments ..
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
type(interval), intent(out) :: interv type(interval), intent(out) :: interv
! .. ! ..
@ -964,12 +984,13 @@ contains
interv%lower = smin interv%lower = smin
interv%upper = smax interv%upper = smax
end subroutine spread_in_coordinate end subroutine spread_in_coordinate
subroutine kdtree2_destroy(tp) subroutine kdtree2_destroy(tp)
! Deallocates all memory for the tree, except input data matrix ! Deallocates all memory for the tree, except input data matrix
! .. Structure Arguments .. ! .. Structure Arguments ..
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
! .. ! ..
call destroy_node(tp%root) call destroy_node(tp%root)
@ -988,6 +1009,7 @@ contains
contains contains
recursive subroutine destroy_node(np) recursive subroutine destroy_node(np)
! .. Structure Arguments .. ! .. Structure Arguments ..
implicit none
type (tree_node), pointer :: np type (tree_node), pointer :: np
! .. ! ..
! .. Intrinsic Functions .. ! .. Intrinsic Functions ..
@ -1007,12 +1029,13 @@ contains
end subroutine destroy_node end subroutine destroy_node
end subroutine kdtree2_destroy end subroutine kdtree2_destroy
subroutine kdtree2_n_nearest(tp,qv,nn,results) subroutine kdtree2_n_nearest(tp,qv,nn,results)
! Find the 'nn' vectors in the tree nearest to 'qv' in euclidean norm ! Find the 'nn' vectors in the tree nearest to 'qv' in euclidean norm
! returning their indexes and distances in 'indexes' and 'distances' ! returning their indexes and distances in 'indexes' and 'distances'
! arrays already allocated passed to this subroutine. ! arrays already allocated passed to this subroutine.
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
real(pReal), target, intent (In) :: qv(:) real(pReal), target, intent (In) :: qv(:)
integer(pInt), intent (In) :: nn integer(pInt), intent (In) :: nn
@ -1050,12 +1073,13 @@ contains
endif endif
! deallocate(sr%pqp) ! deallocate(sr%pqp)
return return
end subroutine kdtree2_n_nearest end subroutine kdtree2_n_nearest
subroutine kdtree2_n_nearest_around_point(tp,idxin,correltime,nn,results) subroutine kdtree2_n_nearest_around_point(tp,idxin,correltime,nn,results)
! Find the 'nn' vectors in the tree nearest to point 'idxin', ! Find the 'nn' vectors in the tree nearest to point 'idxin',
! with correlation window 'correltime', returing results in ! with correlation window 'correltime', returing results in
! results(:), which must be pre-allocated upon entry. ! results(:), which must be pre-allocated upon entry.
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
integer(pInt), intent (In) :: idxin, correltime, nn integer(pInt), intent (In) :: idxin, correltime, nn
type(kdtree2_result), target :: results(:) type(kdtree2_result), target :: results(:)
@ -1093,9 +1117,9 @@ contains
endif endif
deallocate (sr%qv) deallocate (sr%qv)
return return
end subroutine kdtree2_n_nearest_around_point end subroutine kdtree2_n_nearest_around_point
subroutine kdtree2_r_nearest(tp,qv,r2,nfound,nalloc,results) subroutine kdtree2_r_nearest(tp,qv,r2,nfound,nalloc,results)
! find the nearest neighbors to point 'idxin', within SQUARED ! find the nearest neighbors to point 'idxin', within SQUARED
! Euclidean distance 'r2'. Upon ENTRY, nalloc must be the ! Euclidean distance 'r2'. Upon ENTRY, nalloc must be the
! size of memory allocated for results(1:nalloc). Upon ! size of memory allocated for results(1:nalloc). Upon
@ -1106,6 +1130,7 @@ contains
! the smallest ball inside norm r^2 ! the smallest ball inside norm r^2
! !
! Results are NOT sorted unless tree was created with sort option. ! Results are NOT sorted unless tree was created with sort option.
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
real(pReal), target, intent (In) :: qv(:) real(pReal), target, intent (In) :: qv(:)
real(pReal), intent(in) :: r2 real(pReal), intent(in) :: r2
@ -1154,16 +1179,16 @@ contains
endif endif
return return
end subroutine kdtree2_r_nearest end subroutine kdtree2_r_nearest
subroutine kdtree2_r_nearest_around_point(tp,idxin,correltime,r2,& subroutine kdtree2_r_nearest_around_point(tp,idxin,correltime,r2,nfound,nalloc,results)
nfound,nalloc,results)
! !
! Like kdtree2_r_nearest, but around a point 'idxin' already existing ! Like kdtree2_r_nearest, but around a point 'idxin' already existing
! in the data set. ! in the data set.
! !
! Results are NOT sorted unless tree was created with sort option. ! Results are NOT sorted unless tree was created with sort option.
! !
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
integer(pInt), intent (In) :: idxin, correltime, nalloc integer(pInt), intent (In) :: idxin, correltime, nalloc
real(pReal), intent(in) :: r2 real(pReal), intent(in) :: r2
@ -1221,10 +1246,11 @@ contains
deallocate (sr%qv) deallocate (sr%qv)
return return
end subroutine kdtree2_r_nearest_around_point end subroutine kdtree2_r_nearest_around_point
function kdtree2_r_count(tp,qv,r2) result(nfound) function kdtree2_r_count(tp,qv,r2) result(nfound)
! Count the number of neighbors within square distance 'r2'. ! Count the number of neighbors within square distance 'r2'.
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
real(pReal), target, intent (In) :: qv(:) real(pReal), target, intent (In) :: qv(:)
real(pReal), intent(in) :: r2 real(pReal), intent(in) :: r2
@ -1265,13 +1291,13 @@ contains
nfound = sr%nfound nfound = sr%nfound
return return
end function kdtree2_r_count end function kdtree2_r_count
function kdtree2_r_count_around_point(tp,idxin,correltime,r2) & function kdtree2_r_count_around_point(tp,idxin,correltime,r2) result(nfound)
result(nfound)
! Count the number of neighbors within square distance 'r2' around ! Count the number of neighbors within square distance 'r2' around
! point 'idxin' with decorrelation time 'correltime'. ! point 'idxin' with decorrelation time 'correltime'.
! !
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
integer(pInt), intent (In) :: correltime, idxin integer(pInt), intent (In) :: correltime, idxin
real(pReal), intent(in) :: r2 real(pReal), intent(in) :: r2
@ -1315,42 +1341,41 @@ contains
nfound = sr%nfound nfound = sr%nfound
return return
end function kdtree2_r_count_around_point end function kdtree2_r_count_around_point
subroutine validate_query_storage(n) subroutine validate_query_storage(n)
! !
! make sure we have enough storage for n ! make sure we have enough storage for n
! !
use IO, only: IO_error
implicit none
integer(pInt), intent(in) :: n integer(pInt), intent(in) :: n
if (int(size(sr%results,1),pInt) .lt. n) then if (int(size(sr%results,1),pInt) .lt. n) then
write (*,*) 'KD_TREE_TRANS: you did not provide enough storage for results(1:n)' call IO_error (500_pInt,ext_msg='KD_TREE_TRANS: not enough storage for results(1:n)')
stop
return
endif endif
return return
end subroutine validate_query_storage end subroutine validate_query_storage
function square_distance(d, iv,qv) result (res) function square_distance(d, iv,qv) result (res)
! distance between iv[1:n] and qv[1:n] ! distance between iv[1:n] and qv[1:n]
! .. Function Return Value .. ! .. Function Return Value ..
! re-implemented to improve vectorization. ! re-implemented to improve vectorization.
implicit none
real(pReal) :: res real(pReal) :: res
! ..
! ..
! .. Scalar Arguments .. ! .. Scalar Arguments ..
integer(pInt) :: d integer(pInt) :: d
! ..
! .. Array Arguments .. ! .. Array Arguments ..
real(pReal) :: iv(:),qv(:) real(pReal) :: iv(:),qv(:)
! .. ! ..
! .. ! ..
res = sum( (iv(1:d)-qv(1:d))**2 ) res = sum( (iv(1:d)-qv(1:d))**2 )
end function square_distance end function square_distance
recursive subroutine search(node) recursive subroutine search(node)
! !
! This is the innermost core routine of the kd-tree search. Along ! This is the innermost core routine of the kd-tree search. Along
! with "process_terminal_node", it is the performance bottleneck. ! with "process_terminal_node", it is the performance bottleneck.
@ -1358,6 +1383,7 @@ contains
! This version uses a logically complete secondary search of ! This version uses a logically complete secondary search of
! "box in bounds", whether the sear ! "box in bounds", whether the sear
! !
implicit none
type (Tree_node), pointer :: node type (Tree_node), pointer :: node
! .. ! ..
type(tree_node),pointer :: ncloser, nfarther type(tree_node),pointer :: ncloser, nfarther
@ -1424,10 +1450,12 @@ contains
endif endif
endif endif
end if end if
end subroutine search end subroutine search
real(pReal) function dis2_from_bnd(x,amin,amax) result (res) real(pReal) function dis2_from_bnd(x,amin,amax) result (res)
implicit none
real(pReal), intent(in) :: x, amin,amax real(pReal), intent(in) :: x, amin,amax
if (x > amax) then if (x > amax) then
@ -1443,15 +1471,16 @@ contains
endif endif
endif endif
return return
end function dis2_from_bnd end function dis2_from_bnd
logical function box_in_search_range(node, sr) result(res) logical function box_in_search_range(node, sr) result(res)
! !
! Return the distance from 'qv' to the CLOSEST corner of node's ! Return the distance from 'qv' to the CLOSEST corner of node's
! bounding box ! bounding box
! for all coordinates outside the box. Coordinates inside the box ! for all coordinates outside the box. Coordinates inside the box
! contribute nothing to the distance. ! contribute nothing to the distance.
! !
implicit none
type (tree_node), pointer :: node type (tree_node), pointer :: node
type (tree_search_record), pointer :: sr type (tree_search_record), pointer :: sr
@ -1474,14 +1503,15 @@ contains
end do end do
res = .true. res = .true.
return return
end function box_in_search_range end function box_in_search_range
subroutine process_terminal_node(node) subroutine process_terminal_node(node)
! !
! Look for actual near neighbors in 'node', and update ! Look for actual near neighbors in 'node', and update
! the search results on the sr data structure. ! the search results on the sr data structure.
! !
implicit none
type (tree_node), pointer :: node type (tree_node), pointer :: node
! !
real(pReal), pointer :: qv(:) real(pReal), pointer :: qv(:)
@ -1580,14 +1610,15 @@ contains
! !
sr%ballsize = ballsize sr%ballsize = ballsize
end subroutine process_terminal_node end subroutine process_terminal_node
subroutine process_terminal_node_fixedball(node) subroutine process_terminal_node_fixedball(node)
! !
! Look for actual near neighbors in 'node', and update ! Look for actual near neighbors in 'node', and update
! the search results on the sr data structure, i.e. ! the search results on the sr data structure, i.e.
! save all within a fixed ball. ! save all within a fixed ball.
! !
implicit none
type (tree_node), pointer :: node type (tree_node), pointer :: node
! !
real(pReal), pointer :: qv(:) real(pReal), pointer :: qv(:)
@ -1671,13 +1702,14 @@ contains
! Reset sr variables which may have changed during loop ! Reset sr variables which may have changed during loop
! !
sr%nfound = nfound sr%nfound = nfound
end subroutine process_terminal_node_fixedball end subroutine process_terminal_node_fixedball
subroutine kdtree2_n_nearest_brute_force(tp,qv,nn,results) subroutine kdtree2_n_nearest_brute_force(tp,qv,nn,results)
! find the 'n' nearest neighbors to 'qv' by exhaustive search. ! find the 'n' nearest neighbors to 'qv' by exhaustive search.
! only use this subroutine for testing, as it is SLOW! The ! only use this subroutine for testing, as it is SLOW! The
! whole point of a k-d tree is to avoid doing what this subroutine ! whole point of a k-d tree is to avoid doing what this subroutine
! does. ! does.
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
real(pReal), intent (In) :: qv(:) real(pReal), intent (In) :: qv(:)
integer(pInt), intent (In) :: nn integer(pInt), intent (In) :: nn
@ -1710,14 +1742,15 @@ contains
end if end if
end do end do
deallocate (all_distances) deallocate (all_distances)
end subroutine kdtree2_n_nearest_brute_force end subroutine kdtree2_n_nearest_brute_force
subroutine kdtree2_r_nearest_brute_force(tp,qv,r2,nfound,results) subroutine kdtree2_r_nearest_brute_force(tp,qv,r2,nfound,results)
! find the nearest neighbors to 'qv' with distance**2 <= r2 by exhaustive search. ! find the nearest neighbors to 'qv' with distance**2 <= r2 by exhaustive search.
! only use this subroutine for testing, as it is SLOW! The ! only use this subroutine for testing, as it is SLOW! The
! whole point of a k-d tree is to avoid doing what this subroutine ! whole point of a k-d tree is to avoid doing what this subroutine
! does. ! does.
implicit none
type (kdtree2), pointer :: tp type (kdtree2), pointer :: tp
real(pReal), intent (In) :: qv(:) real(pReal), intent (In) :: qv(:)
real(pReal), intent (In) :: r2 real(pReal), intent (In) :: r2
@ -1751,11 +1784,12 @@ contains
call kdtree2_sort_results(nfound,results) call kdtree2_sort_results(nfound,results)
end subroutine kdtree2_r_nearest_brute_force end subroutine kdtree2_r_nearest_brute_force
subroutine kdtree2_sort_results(nfound,results) subroutine kdtree2_sort_results(nfound,results)
! Use after search to sort results(1:nfound) in order of increasing ! Use after search to sort results(1:nfound) in order of increasing
! distance. ! distance.
implicit none
integer(pInt), intent(in) :: nfound integer(pInt), intent(in) :: nfound
type(kdtree2_result), target :: results(:) type(kdtree2_result), target :: results(:)
! !
@ -1767,14 +1801,15 @@ contains
if (nfound .gt. 1_pInt) call heapsort_struct(results,nfound) if (nfound .gt. 1_pInt) call heapsort_struct(results,nfound)
return return
end subroutine kdtree2_sort_results end subroutine kdtree2_sort_results
subroutine heapsort(a,ind,n) subroutine heapsort(a,ind,n)
! !
! Sort a(1:n) in ascending order, permuting ind(1:n) similarly. ! Sort a(1:n) in ascending order, permuting ind(1:n) similarly.
! !
! If ind(k) = k upon input, then it will give a sort index upon output. ! If ind(k) = k upon input, then it will give a sort index upon output.
! !
implicit none
integer(pInt),intent(in) :: n integer(pInt),intent(in) :: n
real(pReal), intent(inout) :: a(:) real(pReal), intent(inout) :: a(:)
integer(pInt), intent(inout) :: ind(:) integer(pInt), intent(inout) :: ind(:)
@ -1826,13 +1861,14 @@ contains
end do end do
a(i)=value; ind(i)=ivalue a(i)=value; ind(i)=ivalue
end do end do
end subroutine heapsort end subroutine heapsort
subroutine heapsort_struct(a,n) subroutine heapsort_struct(a,n)
! !
! Sort a(1:n) in ascending order ! Sort a(1:n) in ascending order
! !
! !
implicit none
integer(pInt),intent(in) :: n integer(pInt),intent(in) :: n
type(kdtree2_result),intent(inout) :: a(:) type(kdtree2_result),intent(inout) :: a(:)
@ -1882,9 +1918,6 @@ contains
end do end do
a(i)=value a(i)=value
end do end do
end subroutine heapsort_struct end subroutine heapsort_struct
end module kdtree2_module end module kdtree2_module
!#############################################################################################################################
! END KDTREE2
!#############################################################################################################################

35
code/prec.f90
View File

@ -19,17 +19,17 @@
!############################################################## !##############################################################
!* $Id$ !* $Id$
!############################################################## !##############################################################
MODULE prec module prec
!############################################################## !##############################################################
implicit none implicit none
private
! *** Precision of real and integer variables *** ! *** Precision of real and integer variables ***
integer, parameter, public :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300 integer, parameter, public :: pReal = selected_real_kind(15,300) ! 15 significant digits, up to 1e+-300
integer, parameter, public :: pInt = selected_int_kind(9) ! up to +- 1e9 integer, parameter, public :: pInt = selected_int_kind(9) ! up to +- 1e9
integer, parameter, public :: pLongInt = 8 ! should be 64bit integer, parameter, public :: pLongInt = 8 ! should be 64bit
real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal
real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-100_pReal real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-100_pReal
! NaN is precision dependent ! NaN is precision dependent
! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html ! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html
@ -43,16 +43,19 @@ real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-100_pReal
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal) real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF0000000000001', pReal)
#endif #endif
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
CONTAINS type, public :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
public :: prec_init
contains
subroutine prec_init subroutine prec_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
implicit none
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
write(6,*) '<<<+- prec init -+>>>' write(6,*) '<<<+- prec init -+>>>'
@ -67,6 +70,6 @@ implicit none
write(6,*) write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
end subroutine end subroutine prec_init
END MODULE prec end module prec

36
code/prec_single.f90
View File

@ -19,17 +19,17 @@
!############################################################## !##############################################################
!* $Id$ !* $Id$
!############################################################## !##############################################################
MODULE prec module prec
!############################################################## !##############################################################
implicit none implicit none
private
! *** Precision of real and integer variables *** ! *** Precision of real and integer variables ***
integer, parameter, public :: pReal = selected_real_kind(6,37) ! 6 significant digits, up to 1e+-37 integer, parameter, public :: pReal = selected_real_kind(6,37) ! 6 significant digits, up to 1e+-37
integer, parameter, public :: pInt = selected_int_kind(9) ! up to +- 1e9 integer, parameter, public :: pInt = selected_int_kind(9) ! up to +- 1e9
integer, parameter, public :: pLongInt = 4 ! should be 64bit integer, parameter, public :: pLongInt = 4 ! should be 64bit
real(pReal), parameter, public :: tol_math_check = 1.0e-5_pReal real(pReal), parameter, public :: tol_math_check = 1.0e-5_pReal
real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-36_pReal real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-36_pReal
! NaN is precision dependent ! NaN is precision dependent
! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html ! from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html
@ -43,15 +43,19 @@ real(pReal), parameter, public :: tol_gravityNodePos = 1.0e-36_pReal
#else #else
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal) real(pReal), parameter, public :: DAMASK_NaN = real(Z'7F800001', pReal)
#endif #endif
type :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
CONTAINS type, public :: p_vec
real(pReal), dimension(:), pointer :: p
end type p_vec
public :: prec_init
contains
subroutine prec_init subroutine prec_init
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
implicit none
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
@ -67,6 +71,6 @@ implicit none
write(6,*) write(6,*)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
end subroutine end subroutine prec_init
END MODULE prec end module prec