Initial commit

CMakeLists.txt

@@ -0,0 +1,51 @@
+cmake_minimum_required(VERSION 2.6 FATAL_ERROR)
+
+set(CMAKE_USER_MAKE_RULES_OVERRIDE ${CMAKE_SOURCE_DIR}/cmake/UserOverride.cmake)
+
+enable_language(Fortran)
+
+project(featom)
+
+set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
+
+set(WITH_LAPACK yes
+    CACHE BOOL "Build with LAPACK support")
+set(WITH_HDF5 no
+    CACHE BOOL "Build with HDF5 support")
+
+# Make sure that CMAKE_BUILD_TYPE is either Debug or Release:
+if (NOT CMAKE_BUILD_TYPE)
+    set(CMAKE_BUILD_TYPE Debug
+        CACHE STRING "Build type (Debug, Release)" FORCE)
+endif ()
+if (NOT (CMAKE_BUILD_TYPE STREQUAL "Debug" OR
+        CMAKE_BUILD_TYPE STREQUAL "Release"))
+    message("${CMAKE_BUILD_TYPE}")
+    message(FATAL_ERROR "CMAKE_BUILD_TYPE must be one of: Debug, Release (current value: '${CMAKE_BUILD_TYPE}')")
+endif ()
+
+if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
+    # gfortran
+    # Enable this if you want to check for single/double corruption (and use
+    # the Debug build):
+    #set(CMAKE_Fortran_FLAGS_DEBUG "${CMAKE_Fortran_FLAGS_DEBUG} -fdefault-real-8")
+endif ()
+
+enable_testing()
+
+add_subdirectory(src)
+add_subdirectory(tests)
+
+message("\n")
+message("Configuration results")
+message("---------------------")
+message("Fortran compiler: ${CMAKE_Fortran_COMPILER}")
+message("Build type: ${CMAKE_BUILD_TYPE}")
+if (CMAKE_BUILD_TYPE STREQUAL "Debug")
+    message("Fortran compiler flags: ${CMAKE_Fortran_FLAGS_DEBUG}")
+else ()
+    message("Fortran compiler flags: ${CMAKE_Fortran_FLAGS_RELEASE}")
+endif ()
+message("Installation prefix: ${CMAKE_INSTALL_PREFIX}")
+message("With LAPACK: ${WITH_LAPACK}")
+message("With HDF5: ${WITH_HDF5}")

LICENSE

@@ -0,0 +1,19 @@
+Copyright (c) 2012 Ondřej Čertík
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.

README.rst

@@ -0,0 +1,32 @@
+Fortran Utilities
+=================
+
+Various Fortran utilities, that can be included into any Fortran
+program.
+
+The modules are mostly independent of each other. Simply copy any modules that
+you need into your project. Tests are in the ``tests`` directory, you can look
+there for examples of usage.
+
+License
+-------
+
+All code is MIT licensed.
+
+Functionality
+-------------
+
+* Types (``dp``)
+* Constants (``pi``, ``e_``, ``i_``)
+* Sorting
+* Saving/loading 2D arrays (``savetxt``, ``loadtxt``)
+* Meshes (exponential, uniform)
+* Cubic splines
+
+Contributors
+------------
+
+List of people who contributed code to this library:
+
+* John E. Pask
+* Ondřej Čertík

cmake/UserOverride.cmake

@@ -0,0 +1,18 @@
+# This overrides the default CMake Debug and Release compiler options.
+# The user can still specify different options by setting the
+# CMAKE_Fortran_FLAGS_[RELEASE,DEBUG] variables (on the command line or in the
+# CMakeList.txt). This files serves as better CMake defaults and should only be
+# modified if the default values are to be changed. Project specific compiler
+# flags should be set in the CMakeList.txt by setting the CMAKE_Fortran_FLAGS_*
+# variables.
+if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
+    # gfortran
+    set(common "-std=f2008 -Wall -Wextra -Wimplicit-interface -fPIC -Werror -fmax-errors=1")
+    set(CMAKE_Fortran_FLAGS_RELEASE_INIT "${common} -O3 -march=native -ffast-math -funroll-loops")
+    set(CMAKE_Fortran_FLAGS_DEBUG_INIT   "${common} -g -fcheck=all -fbacktrace")
+elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
+    # ifort
+    set(common "-std2008 -warn all")
+    set(CMAKE_Fortran_FLAGS_RELEASE_INIT "${common} -xHOST -O3 -no-prec-div -static")
+    set(CMAKE_Fortran_FLAGS_DEBUG_INIT   "${common} -check all")
+endif ()

src/CMakeLists.txt

@@ -0,0 +1,18 @@
+set(SRC
+    constants.f90   sorting.f90           types.f90
+    utils.f90
+    mesh.f90
+)
+
+if(WITH_LAPACK)
+    set(SRC ${SRC} lapack.f90 splines.f90)
+endif()
+
+if(WITH_HDF5)
+    set(SRC ${SRC} h5_utils.f90)
+    # TODO: make this search path more general somehow:
+    # This directory contains the hdf5.mod that we need:
+    include_directories($ENV{SPKG_LOCAL}/include)
+endif()
+
+add_library(fortran_utils ${SRC})

src/constants.f90

@@ -0,0 +1,15 @@
+module constants
+use types, only: dp
+implicit none
+private
+public pi, e_, i_
+
+! Constants contain more digits than double precision, so that
+! they are rounded correctly. Single letter constants contain underscore so
+! that they do not clash with user variables ("e" and "i" are frequently used as
+! loop variables)
+real(dp), parameter :: pi    = 3.1415926535897932384626433832795_dp
+real(dp), parameter :: e_    = 2.7182818284590452353602874713527_dp
+complex(dp), parameter :: i_ = (0, 1)
+
+end module

src/h5_utils.f90

@@ -0,0 +1,124 @@
+module h5_utils
+
+use hdf5, only: HID_T, HSIZE_T, H5T_NATIVE_DOUBLE, H5T_NATIVE_INTEGER, &
+    H5F_ACC_TRUNC_F, h5gcreate_f, h5gclose_f, h5screate_simple_f, &
+    h5sclose_f, h5dcreate_f, h5dclose_f, h5dwrite_f, h5fcreate_f, &
+    h5fclose_f, h5close_f, h5open_f
+use types, only: dp
+use utils, only: stop_error
+implicit none
+private
+public h5_open, h5_close, h5_write_array, h5_file, h5_create_group
+
+type h5_file
+    integer(HID_T) :: file_id
+end type
+
+interface h5_write_array
+    module procedure h5_write_array_int, h5_write_array_real
+end interface
+
+integer :: number_of_open_files = 0
+
+contains
+
+subroutine check(error)
+integer, intent(in) :: error
+if (error /= 0) call stop_error("Error when calling HDF5.")
+end subroutine
+
+type(h5_file) function h5_open(filename)
+character(len=*), intent(in) :: filename
+integer :: error
+integer(HID_T) :: file_id
+if (number_of_open_files == 0) then
+    ! If this is the first call to the HDF5, we need to initialize the
+    ! interface:
+    call h5open_f(error)
+    call check(error)
+end if
+number_of_open_files = number_of_open_files + 1
+call h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, error)
+call check(error)
+h5_open%file_id = file_id
+end function
+
+subroutine h5_close(self)
+type(h5_file), intent(in) :: self
+integer :: error
+call h5fclose_f(self%file_id, error)
+call check(error)
+number_of_open_files = number_of_open_files - 1
+if (number_of_open_files == 0) then
+    ! All files are closed, so we can close the HDF5 interface:
+    call h5close_f(error)
+    call check(error)
+end if
+end subroutine
+
+subroutine h5_write_array_int(self, a_name, A)
+type(h5_file), intent(in) :: self
+character(LEN=*), intent(in) :: a_name
+integer, intent(in) :: A(:)
+
+integer, parameter :: rank = 1
+integer(HID_T) :: dset_id, dspace_id
+integer(HSIZE_T) :: dims(1)
+integer :: error
+
+dims(1) = size(A)
+
+call h5screate_simple_f(rank, dims, dspace_id, error)
+call check(error)
+call h5dcreate_f(self%file_id, a_name, H5T_NATIVE_INTEGER, dspace_id, dset_id, error)
+call check(error)
+
+call h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, A, dims, error)
+call check(error)
+
+call h5dclose_f(dset_id, error)
+call check(error)
+call h5sclose_f(dspace_id, error)
+call check(error)
+end subroutine
+
+subroutine h5_write_array_real(self, a_name, A)
+type(h5_file), intent(in) :: self
+character(LEN=*), intent(in) :: a_name
+real(dp), intent(in) :: A(:)
+
+integer, parameter :: rank = 1
+integer(HID_T) :: dset_id, dspace_id
+integer(HSIZE_T) :: dims(1)
+integer :: error
+
+dims(1) = size(A)
+
+call h5screate_simple_f(rank, dims, dspace_id, error)
+call check(error)
+call h5dcreate_f(self%file_id, a_name, H5T_NATIVE_DOUBLE, dspace_id, dset_id, error)
+call check(error)
+
+call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE, A, dims, error)
+call check(error)
+
+call h5dclose_f(dset_id, error)
+call check(error)
+call h5sclose_f(dspace_id, error)
+call check(error)
+end subroutine
+
+subroutine h5_create_group(self, g_name)
+type(h5_file), intent(in) :: self
+character(LEN=*), intent(in) :: g_name
+
+integer(HID_T) :: group_id
+integer :: error
+
+call h5gcreate_f(self%file_id, g_name, group_id, error)
+call check(error)
+call h5gclose_f(group_id, error)
+call check(error)
+end subroutine
+
+end module

src/lapack.f90

@@ -0,0 +1,100 @@
+module lapack_precision
+implicit none
+private
+public dp
+
+! This is the precision that LAPACK "d" routines were compiled with (typically
+! double precision, unless a special compiler option was used while compiling
+! LAPACK). This "dp" is only used in lapack.f90
+! The "d" routines data type is defined as "double precision", so
+! we make "dp" the same kind as 0.d0 ("double precision"), so
+! as long as LAPACK and this file were compiled with the same compiler options,
+! it will be consistent. (If for example all double precision is promoted to
+! quadruple precision, it will be promoted both in LAPACK and here.)
+integer, parameter:: dp=kind(0.d0)
+
+end module
+
+
+
+
+module lapack
+implicit none
+
+interface
+
+    SUBROUTINE DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
+    use lapack_precision, only: dp
+    INTEGER            INFO, LDA, LDB, N, NRHS
+    INTEGER            IPIV( * )
+    REAL(dp)           A( LDA, * ), B( LDB, * )
+    END SUBROUTINE
+
+    SUBROUTINE DGBSV( N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB, INFO )
+    use lapack_precision, only: dp
+    INTEGER            INFO, KL, KU, LDAB, LDB, N, NRHS
+    INTEGER            IPIV( * )
+    REAL(dp)           AB( LDAB, * ), B( LDB, * )
+    END SUBROUTINE
+
+    SUBROUTINE DSYSV( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, LWORK, INFO )
+    use lapack_precision, only: dp
+    CHARACTER          UPLO
+    INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
+    INTEGER            IPIV( * )
+    REAL(dp)           A( LDA, * ), B( LDB, * ), WORK( * )
+    END SUBROUTINE
+
+    SUBROUTINE DSYSVX( FACT, UPLO, N, NRHS, A, LDA, AF, LDAF, IPIV, B, &
+        LDB, X, LDX, RCOND, FERR, BERR, WORK, LWORK, &
+        IWORK, INFO )
+    use lapack_precision, only: dp
+    CHARACTER          FACT, UPLO
+    INTEGER            INFO, LDA, LDAF, LDB, LDX, LWORK, N, NRHS
+    REAL(dp)           RCOND
+    INTEGER            IPIV( * ), IWORK( * )
+    REAL(dp)           A( LDA, * ), AF( LDAF, * ), B( LDB, * ), &
+                       BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
+    END SUBROUTINE
+
+    SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, &
+                       VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, &
+                       LWORK, IWORK, IFAIL, INFO )
+    use lapack_precision, only: dp
+    CHARACTER          JOBZ, RANGE, UPLO
+    INTEGER            IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
+    REAL(dp)           ABSTOL, VL, VU
+    INTEGER            IFAIL( * ), IWORK( * )
+    REAL(dp)           A( LDA, * ), B( LDB, * ), W( * ), WORK( * ), &
+                       Z( LDZ, * )
+    END SUBROUTINE
+
+    SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, &
+        BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
+    use lapack_precision, only: dp
+    CHARACTER          JOBVL, JOBVR
+    INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
+    REAL(dp)           A( LDA, * ), ALPHAI( * ), ALPHAR( * ), &
+                       B( LDB, * ), BETA( * ), VL( LDVL, * ), &
+                       VR( LDVR, * ), WORK( * )
+    END SUBROUTINE
+
+    SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, &
+        LWORK, IWORK, LIWORK, INFO )
+    use lapack_precision, only: dp
+    CHARACTER          JOBZ, UPLO
+    INTEGER            INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
+    INTEGER            IWORK( * )
+    REAL(dp)           A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
+    END SUBROUTINE
+
+    REAL(dp) FUNCTION DLAMCH( CMACH )
+    use lapack_precision, only: dp
+    CHARACTER          CMACH
+    END FUNCTION
+
+end interface
+
+contains
+
+end module