2

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) ]
--------------------------------------------------^
nicholasr
  • 25
  • 4
  • 1
    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 Feb 23 '17 at 00:26
  • Using explicit kind constants is ugly and non-portable. See http://stackoverflow.com/documentation/fortran/939/data-types/4390/precision-of-floating-point-numbers#t=201702230923398688378 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 Героям слава Feb 23 '17 at 09: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 Feb 23 '17 at 23:19

2 Answers2

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.

jabirali
  • 1,315
  • 7
  • 13
  • 1
    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 http://stackoverflow.com/questions/838310/fortran-90-kind-parameter for a bit of an introduction – Ian Bush Feb 23 '17 at 08:16
  • 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). – Vladimir F Героям слава Feb 23 '17 at 09:18
  • 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.) – jabirali Feb 23 '17 at 18:09
  • 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 Feb 23 '17 at 23:07
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.

nicholasr
  • 25
  • 4
  • 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? – wander95 Mar 01 '17 at 13:14