240304.011511.HKT fortran: introduce ORTHTOL_DFT so that isorth c…

…an be disabled when the compiler is buggy, e.g., NAG Fortran Compiler Release 7.1(Hanzomon) Build 7143 is buggy concerning half-precision floating-point numbers. See Build 7143 is buggy concerning half-precision floating-point numbers

fortran/common/consts.F90

@@ -8,7 +8,7 @@ module consts_mod
 !
 ! Started: July 2020
 !
-! Last Modified: Sunday, February 25, 2024 PM06:00:24
+! Last Modified: Monday, March 04, 2024 AM01:11:04
 !--------------------------------------------------------------------------------------------------!
 
 !--------------------------------------------------------------------------------------------------!
@@ -73,7 +73,7 @@ module consts_mod
 public :: RP, RP_DFT, DP, SP, QP
 public :: ZERO, ONE, TWO, HALF, QUART, TEN, TENTH, PI
 public :: REALMIN, EPS, TINYCV, REALMAX, FUNCMAX, CONSTRMAX, BOUNDMAX
-public :: SYMTOL_DFT
+public :: SYMTOL_DFT, ORTHTOL_DFT
 public :: STDIN, STDOUT, STDERR
 public :: RHOBEG_DFT, RHOEND_DFT, FTARGET_DFT, CTOL_DFT, CWEIGHT_DFT
 public :: ETA1_DFT, ETA2_DFT, GAMMA1_DFT, GAMMA2_DFT
@@ -180,6 +180,16 @@ module consts_mod
 real(RP), parameter :: SYMTOL_DFT = ZERO
 #endif
 
+! ORTHTOL_DFT is the default tolerance for testing orthogonality of matrices.
+! In some cases, due to compiler bugs, we need to disable the test. We signify such cases by setting
+! ORTHTOL_DFT to REALMAX. For instance, NAG Fortran Compiler Release 7.1(Hanzomon) Build 7143 is
+! buggy concerning half-precision floating-point numbers.
+#if (defined __NAG_COMPILER_RELEASE && PRIMA_REAL_PRECISION <= 16) || (PRIMA_RELEASED == 1) || (PRIMA_DEBUGGING == 0)
+real(RP), parameter :: ORTHTOL_DFT = REALMAX
+#else
+real(RP), parameter :: ORTHTOL_DFT = ZERO
+#endif
+
 ! Some default values
 real(RP), parameter :: RHOBEG_DFT = ONE
 real(RP), parameter :: RHOEND_DFT = 1.0E-6_RP

fortran/common/linalg.f90

@@ -39,7 +39,7 @@ module linalg_mod
 !
 ! Started: July 2020
 !
-! Last Modified: Monday, August 14, 2023 PM10:42:59
+! Last Modified: Monday, March 04, 2024 AM01:02:33
 !--------------------------------------------------------------------------------------------------!
 
 implicit none
@@ -1327,7 +1327,7 @@ function isorth(A, tol) result(is_orth)
 !--------------------------------------------------------------------------------------------------!
 ! This function tests whether the matrix A has orthonormal columns up to the tolerance TOL.
 !--------------------------------------------------------------------------------------------------!
-use, non_intrinsic :: consts_mod, only : RP, IK, DEBUGGING
+use, non_intrinsic :: consts_mod, only : RP, IK, DEBUGGING, REALMAX, ORTHTOL_DFT
 use, non_intrinsic :: debug_mod, only : assert
 use, non_intrinsic :: infnan_mod, only : is_nan
 implicit none
@@ -1340,6 +1340,7 @@ function isorth(A, tol) result(is_orth)
 ! Local variables
 character(len=*), parameter :: srname = 'ISORTH'
 integer(IK) :: n
+real(RP) :: tol_loc
 
 ! Preconditions
 if (DEBUGGING) then
@@ -1352,18 +1353,24 @@ function isorth(A, tol) result(is_orth)
 ! Calculation starts !
 !====================!
 
+tol_loc = ORTHTOL_DFT
+if (present(tol)) then
+    tol_loc = tol
+end if
+
 n = int(size(A, 2), kind(n))
 
+! N.B. (20240304): In some cases, due to compiler bugs, we need to disable the test. We signify such
+! cases by setting ORTHTOL_DFT to REALMAX. For instance, NAG Fortran Compiler Release 7.1(Hanzomon)
+! Build 7143 is buggy concerning half-precision floating-point numbers. See the following:
+! https://fortran-lang.discourse.group/t/nagfor-7-1-supports-half-precision-floating-point-numbers-but-with-many-bugs
+is_orth = .true.
 if (n > size(A, 1)) then
     is_orth = .false.
-else if (is_nan(sum(abs(A)))) then
+elseif (is_nan(sum(abs(A)))) then
     is_orth = .false.
-else
-    if (present(tol)) then
-        is_orth = all(abs(matprod(transpose(A), A) - eye(n)) <= max(tol, tol * maxval(abs(A))))
-    else
-        is_orth = all(abs(matprod(transpose(A), A) - eye(n)) <= 0)
-    end if
+elseif (ORTHTOL_DFT < REALMAX) then
+    is_orth = all(abs(matprod(transpose(A), A) - eye(n)) <= max(tol_loc, tol_loc * maxval(abs(A))))
 end if
 
 !====================!
@@ -1844,7 +1851,7 @@ function issymmetric(A, tol) result(is_symmetric)
 is_symmetric = .true.
 if (size(A, 1) /= size(A, 2)) then
     is_symmetric = .false.
-elseif (SYMTOL_DFT < 0.9_RP * REALMAX) then
+elseif (SYMTOL_DFT < REALMAX) then
     is_symmetric = (.not. any(abs(A - transpose(A)) > tol_loc * max(maxval(abs(A)), ONE))) .and. &
         & all(is_nan(A) .eqv. is_nan(transpose(A)))
 end if