Floats causing trouble in Fortran array constructor

Question

I'm a little stuck with some errors related to Fortran array constructors. Maybe you can help me understand the errors? I have two questions in particular (see below).

My minimal working example is:

program debug

use ieee_arithmetic
implicit none

integer :: ii

real, parameter      :: KZERO   = 1
real, parameter      :: LAMBDA  = 1.3

integer, parameter      :: pad     = 5            
integer, parameter      :: nsh     = 60 + 2*pad    
integer, parameter      :: top=nsh-pad, bot=pad+1  

real(kind=8),    parameter :: dt2 = 1.0D-5/2
real(kind=8),    parameter :: VK_SS = 6.0D-10 
real(kind=8),    parameter :: VM_SS = 6.0D-05

real(kind=8), parameter,dimension(nsh) :: k  = [(KZERO*LAMBDA**(ii-pad), ii=1, nsh)] 
real(kind=8), parameter,dimension(nsh) :: NUK_K = [(-dt2*VK_SS*k(ii)**2, ii=1, nsh)]  
real(kind=8), parameter,dimension(nsh) :: NUK_M = [(-dt2*VM_SS*k(ii)**2, ii=1, nsh)]  

real(kind=8),    parameter, dimension(nsh) :: A = [(.0,ii=1,pad), ( REAL(exp(NUK_K(ii))), ii=bot, top), (.0,ii=1,pad)] 
real(kind=8),    parameter, dimension(nsh) :: B = [(.0,ii=1,pad), ( REAL(    NUK_M(ii) ), ii=bot, top), (.0,ii=1,pad)]
!real(kind=8),    parameter, dimension(nsh) :: C = [(.0,ii=1,pad), ( REAL(exp(NUK_M(ii))), ii=bot, top), (.0,ii=1,pad)]

print *,A 
print *,'-------------------------'
print *,B
print *,'-------------------------'
!print *,C

end program debug

First question: How come I need to typecast with REAL() in the constructor where A,B and C are defined? If I do not, I get this error:

/opt/intel/composer_xe_2011_sp1.9.293/bin/intel64/ifort  -o debug debug.f90 
debug.f90(23): error #7113: Each ac-value expression in an array-constructor must have the same type and type parameters.   [EXP]
real(kind=8),    parameter, dimension(nsh) :: A = [(.0,ii=1,pad), ( exp(NUK_K(ii)), ii=bot, top), (.0,ii=1,pad)] 
--------------------------------------------------------------------^

...but the do-loop elements are already of REAL kind=8 since NUK_K is?

Second question: when I uncomment the line defining C (also by an array constructor), I get the following error:

/opt/intel/composer_xe_2011_sp1.9.293/bin/intel64/ifort  -o debug debug.f90 
debug.f90(25): warning #7919: The value was too small when converting to REAL(KIND=4); the result is zero.
real(kind=8),    parameter, dimension(nsh) :: C = [(.0,ii=1,pad), ( REAL(exp(NUK_M(ii))), ii=bot, top), (.0,ii=1,pad)]
--------------------------------------------------------------------^
debug.f90(25): error #7768: This operation on this data type is currently inaccurate.
real(kind=8),    parameter, dimension(nsh) :: C = [(.0,ii=1,pad), ( REAL(exp(NUK_M(ii))), ii=bot, top), (.0,ii=1,pad)]
--------------------------------------------------------------------^

..the warning is OK since I guess it means that the compiler is suggesting exp(-[large number]) -> 0, right? But how do I deal with the error "This operation on this data type is currently inaccurate." ?

I hope you can help me since I've being going at this problem for a long time now!

EDIT

First of all, thank you very much for the helpful answers! I appreciate it a lot. However, the inaccuracy problem remains. Any ideas? My initial guess was to use quads instead of doubles and thereafter cast the result as doubles (quads are not natively supported by my CPU). This does not work either. My new minimal working example is (note I removed the zero-padding to make the example simpler):

program debug

use ieee_arithmetic
use ISO_Fortran_env
implicit none

integer, parameter :: rd = real64
integer, parameter :: rq = real128
integer, parameter :: df = rq ! Default float kind

real(df), parameter      :: KZERO   = 1
real(df), parameter      :: LAMBDA  = 1.3

integer, parameter      :: pad     = 5            
integer, parameter      :: nsh     = 60 + 2*pad    
integer, parameter      :: top=nsh-pad, bot=pad+1  

real(df),    parameter :: dt2 = 1.0D-5/2
real(df),    parameter :: VK_SS = 6.0D-10 
real(df),    parameter :: VM_SS = 6.0D-05

integer :: ii

real(df), parameter,dimension(nsh) :: k  = [(KZERO*LAMBDA**(ii-pad), ii=1, nsh)] 
real(df), parameter,dimension(nsh) :: NUK_K = [(-dt2*VK_SS*k(ii)**2, ii=1, nsh)]  
real(df), parameter,dimension(nsh) :: NUK_M = [(-dt2*VM_SS*k(ii)**2, ii=1, nsh)]  

real(df),    parameter, dimension(nsh) :: A = [ ( exp(NUK_K(ii)) , ii=1, nsh) ] 
real(df),    parameter, dimension(nsh) :: B = [ (     NUK_M(ii)  , ii=1, nsh) ]
real(df),    parameter, dimension(nsh) :: C = [ ( exp(NUK_M(ii)) , ii=1, nsh) ]

print *,A 
print *,'-------------------------'
print *,B
print *,'-------------------------'
print *,C

end program debug

But this still gives the error

/opt/intel/composer_xe_2011_sp1.9.293/bin/intel64/ifort  -o debug debug.f90 
debug.f90(30): error #7768: This operation on this data type is currently inaccurate.
real(df),    parameter, dimension(nsh) :: C = [ ( exp(NUK_M(ii)) , ii=1, nsh) ]
--------------------------------------------------^

The two questions are related, but not in a way that is easy to point you to other resources here. Essentially, think "how are the things on the right evaluated when I ignore what is on the left?". The array constructor doesn't have knowledge of the type of nuk_k; the exponentiation doesn't know about the type of c. So, look up about types of a constructed array, and the result of real(). — francescalus
– francescalus, Commented Feb 23, 2017 at 0:26
Using explicit kind constants is ugly and non-portable. See stackoverflow.com/documentation/fortran/939/data-types/4390/… Furthermore if you have a kind constant say rp, then real(rp) is even shorter then real(kind=8), you can omit the kind=. — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Feb 23, 2017 at 9:24
Thank you for the answers and suggesting the proper way to define different real kinds. I agree, it's better this way. However, the inaccuracy problem remains, which you can see in my edited answer. Any ideas? — nicholasr
– nicholasr, Commented Feb 23, 2017 at 23:19

jabirali · Accepted Answer · 2017-02-23 01:12:10Z

3

First question:

You're trying to concatenate an array (.0,ii=1,pad), which is of the default kind real(4), with the array ( REAL(exp(NUK_K(ii))), ii=bot, top), which is real(8). The easiest solution is to make the first array into a real(8) array too, by e.g. replacing the floating-point number .0 with 0.0d0 (scientific notation for double-precision numbers).

Example:

real(kind=8), parameter, dimension(nsh) :: A = [(0.0d0,ii=1,pad), ( exp(NUK_K(ii)), ii=bot, top), (0.0d0,ii=1,pad)]
real(kind=8), parameter, dimension(nsh) :: B = [(0.0d0,ii=1,pad), (     NUK_M(ii) , ii=bot, top), (0.0d0,ii=1,pad)]
real(kind=8), parameter, dimension(nsh) :: C = [(0.0d0,ii=1,pad), ( exp(NUK_M(ii)), ii=bot, top), (0.0d0,ii=1,pad)]

Second question:

By default, real(num) doesn't typecast the number num to real(8), but to the default kind real(4). If you want to typecast it to real(8), you have to do so explicitly, using the notation real(num,kind=8).

So in other words, what you have been doing, is explicitly typecasting ( exp(NUK_M(ii)), ii=bot, top) from double-precision to single-precision float, then concatenating it with a single-precision array of zeros, and then typecasting the result to a double-precision array again. That's why the compiler was complaining about loss of precision.

answered Feb 23, 2017 at 1:12

jabirali

1,3467 silver badges13 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

Ian Bush Over a year ago

You have no guarantee that 0.0d0 is a real with kind value 8, and I can show you cases where this is not true. Further you have no guarantee that a kind vale of 8 for reals is even valid for the compiler you are using, again I can show an example. In this day and age you shouldn't be using d0, and you should never use explicit constants for kind values. Rather see stackoverflow.com/questions/838310/fortran-90-kind-parameter for a bit of an introduction

Vladimir F Героям слава Over a year ago

Default kind is a default kind, only in some compilers it is 4. It is 1 in several other compilers which you can encounter in the real world. Furthermore, most compilers can be set to use 8-byte reals as the default kind and then double precision (or 1d0) will be 16-byte although kind=8 will still be 8-byte in gfortran. In other words gfortran can be set to use kind=8 as double precision and some people do use this setting (even in questions here).

jabirali Over a year ago

The reason I did not use kind-variables in this answer, was that I was trying to address the OPs problems as specifically as possible given the example code he/she provided. But I agree that in general, using kind-variables is a better practice. (In my own code, I do always define dp = real64 where real64 is imported from iso_fortran_env, and would write 0.0_dp instead of 0.0d0 for the constants.)

nicholasr Over a year ago

Thank you very much for the answers. They were very helpful, and I have adopted the shorthand approach for defining float kinds. However, the inaccuracy problem remains, as you can see from me edited question above. Any ideas?

nicholasr · Accepted Answer · 2017-03-01 13:07:00Z

0

Followup for those who might be interested:

I decided to do a hacky solution to address the problem regarding the precision of the intrinsic function exp(). The problem is that I wish to calculate exp(-[very large number]), which becomes too close to zero for even quad floats to represent the large (negative) exponent. The problem is, that I can not simply set these entries to zero since that messes with the numeric of my full problem. What I therefore instead have done is to replace exp(-[number lager than 7.0D2]) with exp(-7.0D2), which is close to the limit of what is representable using doubles.

answered Mar 1, 2017 at 13:07

nicholasr

254 bronze badges

1 Comment

wander95 Over a year ago

Why not define your own exp_new(x) that calls intrinsic exp(x) for x< 7.0D2 and uses a simple Taylor expansion for x> 7.0D2?

Collectives™ on Stack Overflow

Floats causing trouble in Fortran array constructor

EDIT

2 Answers 2

4 Comments

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

EDIT

2 Answers 2

4 Comments

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related