Default REAL and DOUBLE PRECISION - Results in Incorrectness of Results

[brianlamm]

I will most likely have to file a problem report, but want to see if anyone else has experienced this: where building in "debug", I set "Default Real KIND" to 8 and "Default Double Precision KIND" to 16. I link to Intel MKL as well.

The results of the application built with above settings yields significantly incorrect results. Changing those two above settings back to "" (which sets the above to 4 and 8, respectively), the built application delivers as expected.

Those builds are for x64, host machine is x64 Windows.

I have not yet experimented with "release" for x64, and have not experimented with "release" or "debug" building for ia32.

I know that's high level (read vague), but I'm interested in finding out if anyone else has experienced that resulting incorrectness of results.

Brian Lamm

[Steven_L_Intel1]

You say you link to Intel MKL - how do you select the proper MKL routine to call for the data kind?

[brianlamm]

Oh, that's absolutely simple, robust, and portable across any f95-compliant compiler. I use f95 interfaces to MKL, which call out the real kind as WP = KIND(1.0D0) for double precison routines (WP is declared as a parameter).

In my code, I do the same simple, robust, completely portable implementation: I declare a parameter DP = KIND(1.0D0) for the kind type parameter for double precision reals.

Therefore the default KIND is taken from whatever default is in place for single and double precision reals at the time of the build.

Nowhere in any of my or any of MKL interfaces are explicit numerical values referenced, it's always WP = KIND(1.0E0) for single precision and WP = KIND(1.0D0) for double precision, so that the default kind numbers in affect at the time of build are used.

I do what Intel MKL documentation calls what an "advanced user" does: I do not use MKL capability to pre-build f95 interfaces, I do it directly by INCLUDE'ing the correct source files. So when IVF builder gets to the MKL definitions, they are included as source and the "WP" parameters get set to what IVF has defined as "default".

Simple.

Brian Lamm

[brianlamm]

Here's the relevant code from my app:

!########################################################################################################
!
MODULE mod_ircl_intrinsic_types
PUBLIC
INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9)
INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4)
INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2)
INTEGER, PARAMETER :: SP = KIND(1.0)
INTEGER, PARAMETER :: DP = KIND(1.0D0)
INTEGER, PARAMETER :: SPC = KIND((1.0,1.0))
INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0))
INTEGER, PARAMETER :: LGL = KIND(.TRUE.)
REAL(SP), PARAMETER :: PI = 3.141592653589793238462643383279502884197_sp
REAL(SP), PARAMETER :: PIO2 = 1.57079632679489661923132169163975144209858_sp
REAL(SP), PARAMETER :: TWOPI = 6.283185307179586476925286766559005768394_sp
REAL(DP), PARAMETER :: PI_D = 3.141592653589793238462643383279502884197_dp
REAL(DP), PARAMETER :: PIO2_D = 1.57079632679489661923132169163975144209858_dp
REAL(DP), PARAMETER :: TWOPI_D = 6.283185307179586476925286766559005768394_dp
END MODULE mod_ircl_intrinsic_types
!
!########################################################################################################
!
! The following is Intel-MKL-specific: the file INCLUDEd below houses in a MODULE
! the INTEFACEs of the Fortran-77 "versions" of Intel-MKL, and are USEd by the
! Fortran-95 interface versions, the source codes of which are INCLUDEd in
! MODULE mod_BLAS_realcplx_int.
!
! NOTE: The following INCLUDE must be placed before MODULE mod_BLAS_realcmplx_int
! since the MODULE PROCEDUREs called out in MODULE mod_BLAS_realcmplx_int
! in the generic INTERFACEs have source files (in the INCLUDEs) that reference
! the MODULE MKL77_BLAS found in 'mkl_blas_protos.f90'
!
INCLUDE 'mkl_blas_protos.f90' !* ( one MODULE (interfaces only) : to date (100406) interfaces
!* do not explicitly appear in Lotek -code )
!
!
!########################################################################################################
!
MODULE mod_BLAS_realcplx_int
!
PUBLIC
!
! import BLAS through Intel MKL ( MKL-specific )
! NOTE: updated to MKL Fortran 95 interfaces 12/05/2005
!
INTERFACE DOT
MODULE PROCEDURE SDOT_MKL95, DDOT_MKL95
END INTERFACE DOT
!
INTERFACE GEMV
MODULE PROCEDURE SGEMV_MKL95, DGEMV_MKL95, &
CGEMV_MKL95, ZGEMV_MKL95
END INTERFACE GEMV

INTERFACE NRM2
MODULE PROCEDURE SNRM2_MKL95, DNRM2_MKL95, &
SCNRM2_MKL95, DZNRM2_MKL95
END INTERFACE NRM2
!
CONTAINS
!
INCLUDE 'sdot.f90'
INCLUDE 'ddot.f90'
!
INCLUDE 'sgemv.f90'
INCLUDE 'dgemv.f90'
INCLUDE 'cgemv.f90'
INCLUDE 'zgemv.f90'
!
INCLUDE 'snrm2.f90'
INCLUDE 'dnrm2.f90'
INCLUDE 'scnrm2.f90'
INCLUDE 'dznrm2.f90'
!
!-------------------------------------------------------------------------------------------------------
!
END MODULE mod_BLAS_realcplx_int
!
!########################################################################################################

so the correct real kinds are used because I have defined generics which disambiguate based on KIND alone.

A snippet of the module "MODULE MKL77_BLAS" (MKL-provided) housed in file "mkl_blas_protos.f90" (MKL-provided) for the "dot" function is

MODULE MKL77_BLAS
.
.
INTERFACE MKL77_DOT
PURE FUNCTION SDOT(N,X,INCX,Y,INCY)
INTEGER, PARAMETER :: WP = KIND(1.0E0)
REAL(WP) :: SDOT
REAL(WP), INTENT(IN) :: X(*)
REAL(WP), INTENT(IN) :: Y(*)
INTEGER, INTENT(IN) :: INCX
INTEGER, INTENT(IN) :: INCY
INTEGER, INTENT(IN) :: N
END FUNCTION SDOT
PURE FUNCTION DDOT(N,X,INCX,Y,INCY)
INTEGER, PARAMETER :: WP = KIND(1.0D0)
REAL(WP) :: DDOT
REAL(WP), INTENT(IN) :: X(*)
REAL(WP), INTENT(IN) :: Y(*)
INTEGER, INTENT(IN) :: INCX
INTEGER, INTENT(IN) :: INCY
INTEGER, INTENT( IN) :: N
END FUNCTION DDOT
END INTERFACE MKL77_DOT
.
.

The source code in sdot.f90 (MKL-provided) is as

PURE FUNCTION SDOT_MKL95(X,Y)
! MKL Fortran77 call:
! SDOT(N,X,INCX,Y,INCY)
! <<< Use statements >>>
USE MKL77_BLAS, ONLY: MKL77_DOT
! <<< Implicit statement >>>
IMPLICIT NONE
! <<< Kind parameter >>>
INTEGER, PARAMETER :: WP = KIND(1.0E0)
REAL(WP) :: SDOT_MKL95
! <<< Array arguments >>>
REAL(WP), INTENT(IN) :: X(:)
REAL(WP), INTENT(IN) :: Y(:)
! <<< Local declarations >>>
! <<< Parameters >>>
CHARACTER(LEN=3), PARAMETER :: SRNAME = 'DOT'
! <<< Local scalars >>>
INTEGER :: INCX
INTEGER :: INCY
INTEGER :: N
! <<< Intrinsic functions >>>
INTRINSIC SIZE
! <<< Executable statements >>>
! <<< Init optional and skipped scalars >>>
INCX = 1
INCY = 1
N = SIZE(X)
! <<< Call blas77 routine >>>
SDOT_MKL95 = MKL77_DOT(N,X,INCX,Y,INCY)
END FUNCTION SDOT_MKL95

and the source code in "ddot.f90" (MKL-provided) is as

PURE FUNCTION DDOT_MKL95(X,Y)
! MKL Fortran77 call:
! DDOT(N,X,INCX,Y,INCY)
! <<< Use statements >>>
USE MKL77_BLAS, ONLY: MKL77_DOT
! <<< Implicit statement >>>
IMPLICIT NONE
! <<< Kind parameter >>>
INTEGER, PARAMETER :: WP = KIND(1.0D0)
REAL(WP) :: DDOT_MKL95
! <<< Array arguments >>>
REAL(WP), INTENT(IN) :: X(:)
REAL(WP), INTENT(IN) :: Y(:)
! <<< Local declarations >>>
! <<< Parameters >>>
CHARACTER(LEN=3), PARAMETER :: SRNAME = 'DOT'
! <<< Local scalars >>>
INTEGER :: INCX
INTEGER :: INCY
INTEGER :: N
! <<< Intrinsic functions >>>
INTRINSIC SIZE
! <<< Executable statements >>>
! <<< Init optional and skipped scalars >>>
INCX = 1
INCY = 1
N = SIZE(X)
! <<< Call blas77 routine >>>
DDOT_MKL95 = MKL77_DOT(N,X,INCX,Y,INCY)
END FUNCTION DDOT_MKL95

so, as you can see, the "working precision" is set at *build* time, which, as I wrote above, where IVF encounters those definitions, IVF ha s been set to make the default kind type parameters for real single and double precision to whatever I "told" IVF to set them to, in the "incorrect results scenario", they have been respectively set to 8 and 16.

When I "tell" IVF to set default single and double precision reals to 4 and 8, expected (correct) results are obtained.

Brian Lamm

edit: and, besides, I do *not* get the incorrect results when I set IVF to make default INTEGER kind 8, and there sure are integers used in the MKL routines I reference.

[Steven_L_Intel1]

Thanks for the detail. I suggest you find out where in the computations the results start to diverge.

[brianlamm]

Well, the verdict is in. MKL has no support for 16-byte anything. It cannot support single precision defaulting to 64-bit, and it cannot support double precision defaulting to 128-bit. MKL can accommodate 64-bit integers, but the ILP64 version of MKL would have to be used.

That's per MKL moderator Tim.

Seems like a needless roadblock, however. Why not compile "double" versions (64-bit single and 128-bit double precision) and make them available?

Brian Lamm

[Steven_L_Intel1]

I suspect that part of the problem is that MKL's history is as a performance library and 128-bit floating is emulated in software AND is Fortran-only (so far in our implementations). It would be largely pointless to try to optimize REAL(16) versions of MKL routines, though I can understand that having something available at all would be a convenience.

Please file a feature request with Intel Premier Support (MKL product) requesting REAL(16) versions of MKL routines and explain your reasons for wanting it. In the meantime, you might do well by finding Netlib versions in Fortran and compiling them for REAL(16).

[brianlamm]

Hi Steve,

"largely pointless" ... ? Why would it be pointless?

Thanks for the advice on compiling Netlib versions for REAL(16).

Unfortunately, however, MKL pervades the application I'm developing,apart of which I exposed above, and I use Pardiso and the "dss" routines.

I'll check to see if there's a PARDISO available in source. But even if I could find such, I am somewhat dubious I and IVF would do as excellent a job of optimization as what's found in MKL.

Brian Lamm

[Steven_L_Intel1]

I say it's pointless because any optimization done in the algorithm by the MKL programmers would be swamped by the time taken to do the REAL(16) operations in the software library. It would also mean having to write special versions of the MKL routines that look very different from the single/double versions, with limited payback and increased support cost.

Do you really need REAL(16) for this part of the application? Performance of REAL(16) operations can be 20-50 times slower than double-precision.

[brianlamm]

Yeah, there's one "critical" area accessed every "n" iterations (configurable, starts as user input and then later determines for itself) in which the "quad" precision would be useful for error esitmates on calculations carried out in "double" precision.

I'll find my way around this (need/desire for precise error estimates of "double" precision calcs).

Thanks for the clarifications.

Brian Lamm