Merge branch '460-suppression-of-division-by-zero-when-computing-eige…

…nvalues-in-math_eigh33' into 'development'

avoid division by zero

Closes #460

See merge request damask/DAMASK!992

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,23 @@ 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 = 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(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
 
 end module math

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)::]