From 4583c17080e1ce956dc5778632b1ab09d7c22436 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Dec 2021 07:00:28 +0100 Subject: [PATCH 1/5] corrent 'intent' specification - http://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga5ee879032a8365897c3ba91e3dc8d512.html - http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga56d9c860ce4ce42ded7f914fdb0683ff.html --- src/LAPACK_interface.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90 index 7d3043ed0..91d2451a4 100644 --- a/src/LAPACK_interface.f90 +++ b/src/LAPACK_interface.f90 @@ -23,7 +23,7 @@ module LAPACK_interface integer, intent(in) :: n,nrhs,lda,ldb real(pReal), intent(inout), dimension(lda,n) :: a integer, intent(out), dimension(n) :: ipiv - real(pReal), intent(out), dimension(ldb,nrhs) :: b + real(pReal), intent(inout), dimension(ldb,nrhs) :: b integer, intent(out) :: info end subroutine dgesv @@ -39,7 +39,7 @@ module LAPACK_interface use prec integer, intent(in) :: n,lda,lwork real(pReal), intent(inout), dimension(lda,n) :: a - integer, intent(out), dimension(n) :: ipiv + integer, intent(in), dimension(n) :: ipiv real(pReal), intent(out), dimension(max(1,lwork)) :: work integer, intent(out) :: info end subroutine dgetri From 59bb264b5fd9d3d7dc7d87b8aedba66727506d78 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Dec 2021 07:09:52 +0100 Subject: [PATCH 2/5] LAPACK routines can be considered pure all arguments have 'intent' specification and don't access any global variables. output to screen only occurs in the case that someting goes wrong --- src/LAPACK_interface.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/LAPACK_interface.f90 b/src/LAPACK_interface.f90 index 91d2451a4..e11f74875 100644 --- a/src/LAPACK_interface.f90 +++ b/src/LAPACK_interface.f90 @@ -6,7 +6,7 @@ module LAPACK_interface interface - subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info) + pure subroutine dgeev(jobvl,jobvr,n,a,lda,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info) use prec character, intent(in) :: jobvl,jobvr integer, intent(in) :: n,lda,ldvl,ldvr,lwork @@ -18,7 +18,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgeev - subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info) + pure subroutine dgesv(n,nrhs,a,lda,ipiv,b,ldb,info) use prec integer, intent(in) :: n,nrhs,lda,ldb real(pReal), intent(inout), dimension(lda,n) :: a @@ -27,7 +27,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgesv - subroutine dgetrf(m,n,a,lda,ipiv,info) + pure subroutine dgetrf(m,n,a,lda,ipiv,info) use prec integer, intent(in) :: m,n,lda real(pReal), intent(inout), dimension(lda,n) :: a @@ -35,7 +35,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgetrf - subroutine dgetri(n,a,lda,ipiv,work,lwork,info) + pure subroutine dgetri(n,a,lda,ipiv,work,lwork,info) use prec integer, intent(in) :: n,lda,lwork real(pReal), intent(inout), dimension(lda,n) :: a @@ -44,7 +44,7 @@ module LAPACK_interface integer, intent(out) :: info end subroutine dgetri - subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info) + pure subroutine dsyev(jobz,uplo,n,a,lda,w,work,lwork,info) use prec character, intent(in) :: jobz,uplo integer, intent(in) :: n,lda,lwork From fb51e3c4cd45372b394f5963c88063472a3a776a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 29 Dec 2021 07:19:26 +0100 Subject: [PATCH 3/5] functions have no side-effects, hence 'pure' --- src/lattice.f90 | 4 ++-- src/math.f90 | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index 50d4b6c51..3089aa30b 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -2070,7 +2070,7 @@ end function getlabels !> @brief Equivalent Poisson's ratio (ν) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- -function lattice_equivalent_nu(C,assumption) result(nu) +pure function lattice_equivalent_nu(C,assumption) result(nu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) @@ -2103,7 +2103,7 @@ end function lattice_equivalent_nu !> @brief Equivalent shear modulus (μ) !> @details https://doi.org/10.1143/JPSJ.20.635 !-------------------------------------------------------------------------------------------------- -function lattice_equivalent_mu(C,assumption) result(mu) +pure function lattice_equivalent_mu(C,assumption) result(mu) real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation) character(len=5), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress) diff --git a/src/math.f90 b/src/math.f90 index 2d8f8fb3b..34cb6eecb 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -512,7 +512,7 @@ end subroutine math_invert33 !-------------------------------------------------------------------------------------------------- !> @brief Inversion of symmetriced 3x3x3x3 matrix !-------------------------------------------------------------------------------------------------- -function math_invSym3333(A) +pure function math_invSym3333(A) real(pReal),dimension(3,3,3,3) :: math_invSym3333 @@ -538,7 +538,7 @@ end function math_invSym3333 !-------------------------------------------------------------------------------------------------- !> @brief invert quadratic matrix of arbitrary dimension !-------------------------------------------------------------------------------------------------- -subroutine math_invert(InvA, error, A) +pure subroutine math_invert(InvA, error, A) real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA @@ -996,7 +996,7 @@ end subroutine math_normal !-------------------------------------------------------------------------------------------------- !> @brief eigenvalues and eigenvectors of symmetric matrix !-------------------------------------------------------------------------------------------------- -subroutine math_eigh(w,v,error,m) +pure subroutine math_eigh(w,v,error,m) real(pReal), dimension(:,:), intent(in) :: m !< quadratic matrix to compute eigenvectors and values of real(pReal), dimension(size(m,1)), intent(out) :: w !< eigenvalues @@ -1021,7 +1021,7 @@ end subroutine math_eigh !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @details See http://arxiv.org/abs/physics/0610206 (DSYEVH3) !-------------------------------------------------------------------------------------------------- -subroutine math_eigh33(w,v,m) +pure subroutine math_eigh33(w,v,m) real(pReal), dimension(3,3),intent(in) :: m !< 3x3 matrix to compute eigenvectors and values of real(pReal), dimension(3), intent(out) :: w !< eigenvalues @@ -1114,7 +1114,7 @@ end function math_rotationalPart !> @brief Eigenvalues of symmetric matrix ! will return NaN on error !-------------------------------------------------------------------------------------------------- -function math_eigvalsh(m) +pure function math_eigvalsh(m) real(pReal), dimension(:,:), intent(in) :: m !< symmetric matrix to compute eigenvalues of real(pReal), dimension(size(m,1)) :: math_eigvalsh @@ -1137,7 +1137,7 @@ end function math_eigvalsh !> but apparently more stable solution and has general LAPACK powered version for arbritrary sized !> matrices as fallback !-------------------------------------------------------------------------------------------------- -function math_eigvalsh33(m) +pure function math_eigvalsh33(m) real(pReal), intent(in), dimension(3,3) :: m !< 3x3 symmetric matrix to compute eigenvalues of real(pReal), dimension(3) :: math_eigvalsh33,I From 2f74e0d0700b1f581226c17daceeb4d2a3082f42 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 Dec 2021 21:53:00 +0100 Subject: [PATCH 4/5] avoid failing self test increase number of samples to have less corner cases. Needs to be allocatable to avoid stack/heap issue on ifort --- src/math.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 34cb6eecb..3ea4d6a08 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1432,9 +1432,11 @@ subroutine selfTest error stop 'math_LeviCivita' normal_distribution: block - real(pReal), dimension(500000) :: r + integer, parameter :: N = 1000000 + real(pReal), dimension(:), allocatable :: r real(pReal) :: mu, sigma + allocate(r(N)) call random_number(mu) call random_number(sigma) @@ -1443,11 +1445,11 @@ subroutine selfTest call math_normal(r,mu,sigma) - if (abs(mu -sum(r)/real(size(r),pReal))>5.0e-2_pReal) & + if (abs(mu -sum(r)/real(N,pReal))>5.0e-2_pReal) & error stop 'math_normal(mu)' - mu = sum(r)/real(size(r),pReal) - if (abs(sigma**2 -1.0_pReal/real(size(r)-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) & + mu = sum(r)/real(N,pReal) + if (abs(sigma**2 -1.0_pReal/real(N-1,pReal) * sum((r-mu)**2))/sigma > 5.0e-2_pReal) & error stop 'math_normal(sigma)' end block normal_distribution From b34655b7fc3cc2bc999af482a6e1abd828d7e812 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 30 Dec 2021 21:55:49 +0100 Subject: [PATCH 5/5] functions without side-effects are 'pure' basically all 'getter' functions should be pure --- src/phase_mechanical.f90 | 6 +++--- src/phase_mechanical_elastic.f90 | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 37e6fe7bf..1917d81c9 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -168,17 +168,17 @@ submodule(phase) mechanical integer, intent(in) :: ph,en end function plastic_dislotwin_homogenizedC - module function elastic_C66(ph,en) result(C66) + pure module function elastic_C66(ph,en) result(C66) real(pReal), dimension(6,6) :: C66 integer, intent(in) :: ph, en end function elastic_C66 - module function elastic_mu(ph,en) result(mu) + pure module function elastic_mu(ph,en) result(mu) real(pReal) :: mu integer, intent(in) :: ph, en end function elastic_mu - module function elastic_nu(ph,en) result(nu) + pure module function elastic_nu(ph,en) result(nu) real(pReal) :: nu integer, intent(in) :: ph, en end function elastic_nu diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index acbbac236..24797a880 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -86,7 +86,7 @@ end subroutine elastic_init !-------------------------------------------------------------------------------------------------- !> @brief return 6x6 elasticity tensor !-------------------------------------------------------------------------------------------------- -module function elastic_C66(ph,en) result(C66) +pure module function elastic_C66(ph,en) result(C66) integer, intent(in) :: & ph, & @@ -140,7 +140,7 @@ end function elastic_C66 !-------------------------------------------------------------------------------------------------- !> @brief return shear modulus !-------------------------------------------------------------------------------------------------- -module function elastic_mu(ph,en) result(mu) +pure module function elastic_mu(ph,en) result(mu) integer, intent(in) :: & ph, & @@ -157,7 +157,7 @@ end function elastic_mu !-------------------------------------------------------------------------------------------------- !> @brief return Poisson ratio !-------------------------------------------------------------------------------------------------- -module function elastic_nu(ph,en) result(nu) +pure module function elastic_nu(ph,en) result(nu) integer, intent(in) :: & ph, &