From 066c494c4bae8b6958eabad8dcb055877d28802d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 23 Oct 2024 23:26:37 +0200 Subject: [PATCH 1/2] avoid division by zero happens for all-zero tensor (e.g. during initial call to dislotwin/shear banding). --- src/math.f90 | 16 +++++++++++++++- src/prec.f90 | 2 +- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index c369eeea8..1e48d8cd7 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1038,7 +1038,7 @@ pure subroutine math_eigh33(w,v,m) T = maxval(abs(w)) U = max(T, T**2) - threshold = sqrt(5.68e-14_pREAL * U**2) + threshold = max(sqrt(5.68e-14_pREAL * U**2),PREAL_MIN) v(1:3,1) = [m(1,3)*w(1) + v(1,2), & m(2,3)*w(1) + v(2,2), & @@ -1468,6 +1468,20 @@ subroutine math_selfTest() error stop 'math_normal(sigma)' end block normal_distribution + t33 = 0.0_pREAL + call math_eigh33(v3_1,t33_2,t33) + if (any(dNeq0(v3_1))) error stop 'math_eigh33/zero eigenvalues (values)' + if (any(dNeq(t33_2,math_I3))) error stop 'math_eigh33/zero eigenvalues (vectors)' + t33 = math_I3 + call random_number(r) + d = int(r*2.0_pREAL) + 1 + t33(d,d) = 5.0_pREAL + r*10.0_pREAL + t33(mod(d,3)+1,mod(d,3)+1) = 20.0_pREAL + r*10.0_pREAL + call math_eigh33(v3_1,t33_2,t33) + if (any(dNeq(v3_1,[1.0_pREAL,t33(d,d),t33(mod(d,3)+1,mod(d,3)+1)]))) & + error stop 'math_eigh33/non-zero eigenvalues (values)' + if (any(dNeq(t33_2,math_I3))) error stop 'math_eigh33/non-zero eigenvalues (vectors)' + end subroutine math_selfTest end module math diff --git a/src/prec.f90 b/src/prec.f90 index f9cb6052d..610465707 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -35,7 +35,7 @@ module prec real(pREAL), private, parameter :: PREAL_EPSILON = epsilon(0.0_pREAL) !< minimum positive number such that 1.0 + EPSILON /= 1.0. - real(pREAL), private, parameter :: PREAL_MIN = tiny(0.0_pREAL) !< smallest normalized floating point number + real(pREAL), public, parameter :: PREAL_MIN = tiny(0.0_pREAL) !< smallest normalized floating point number integer, dimension(0), parameter :: emptyIntArray = [integer::] real(pREAL), dimension(0), parameter :: emptyRealArray = [real(pREAL)::] From e820e3454d41602c5e935d4c6feeda1498dc54b5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Oct 2024 08:44:24 +0200 Subject: [PATCH 2/2] consider ordering of eigenvectors --- src/math.f90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index 1e48d8cd7..3156afef9 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1472,15 +1472,18 @@ subroutine math_selfTest() call math_eigh33(v3_1,t33_2,t33) if (any(dNeq0(v3_1))) error stop 'math_eigh33/zero eigenvalues (values)' if (any(dNeq(t33_2,math_I3))) error stop 'math_eigh33/zero eigenvalues (vectors)' + t33 = math_I3 call random_number(r) - d = int(r*2.0_pREAL) + 1 + d = nint(r*2.0_pREAL) + 1 t33(d,d) = 5.0_pREAL + r*10.0_pREAL t33(mod(d,3)+1,mod(d,3)+1) = 20.0_pREAL + r*10.0_pREAL call math_eigh33(v3_1,t33_2,t33) if (any(dNeq(v3_1,[1.0_pREAL,t33(d,d),t33(mod(d,3)+1,mod(d,3)+1)]))) & error stop 'math_eigh33/non-zero eigenvalues (values)' - if (any(dNeq(t33_2,math_I3))) error stop 'math_eigh33/non-zero eigenvalues (vectors)' + if (any(dNeq(math_I3(1:3,mod(d+1,3)+1),t33_2(1:3,1)))) error stop 'math_eigh33/min eigenvector' + if (any(dNeq(math_I3(1:3,d ),t33_2(1:3,2)))) error stop 'math_eigh33/mid eigenvector' + if (any(dNeq(math_I3(1:3,mod(d,3)+1 ),t33_2(1:3,3)))) error stop 'math_eigh33/max eigenvector' end subroutine math_selfTest