- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

I've been having some trouble with choosing the variables precision in my code...

For a while i've been declaring variables as real(kind=8) :: var and I now understand that this is not very portable and has some other complications (which I don't really understand fully) but basically I'm getting a lot of imprecision in numerical calculations. Now I was thinking of adding a module where I would set:

INTEGER, PARAMETER :: R8 = SELECTED_REAL_KIND (15, 300)

And variable declaration should now be real(R8) :: var1,var2.

Here follows the questions:

- I usually initialized variables as var1 = 1.0D0 and now I'm using var1 = 1.0_R8... What should I do with var1 = 1.0D-20?

- Doing numerical calculations, I've found that real(var1 - var2) will give me an imprecise result and using real(var1-var2,kind=R8) gives me a more precise result (I know that using "kind=" is redundant). Is there a better way to increase the precision of numerical calculations, even by enabling any flag in the compiler? I already use -mp1 -fpmodel strict and -prec-div

- Normally I define complex variables as cmplx(0.0_R8,0.0_R8,kind=R8). I believe that "kind=R8" is redundant in this case but I really want to ensure that the result is in double precision. Should I use dcmplx?

- Final question, I also use fftw3 library which uses some constants from iso_c_binding environment and I have to convert some variables to the appropriate kind (ex: from R8 to C_DOUBLE). Should I always use the types from iso_c_binding or should I stick with some minor conversions in the code?

Thank you for your help!

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

First bullet point: var1=1.0E-20_R8 . You have to change the D to an E as well as appending the _R8. The decision could have been made that the kind overrode the D exponent, but it wasn't.

Second bullet point: the kind=R8 is redundant only if at least one of var1 or var2 is of type COMPLEX and at least one has kind=R8. If var1-var2 is not of type COMPLEX then REAL(var1-var2) will be default REAL: kind(REAL((var1-var2)) = kind(1.0) which probably is not what you want. If var1-var2 is of type COMPLEX, then REAL(var1-var2) just takes the real part of var1-var2, leaving its KIND intact: kind(REAL(var1-var2)) = kind(var1-var2).

Third bullet point: this is the most counterintuitive case I think. You pretty much always need the kind=R8 here, otherwise the result will be default COMPLEX: kind(cmplx(0.0_R8,0.0_R8,kind=R8)) = R8, but kind(cmplx(0.0_R8,0.0_R8)) = kind(1.0). Don't use dcmplx because it makes it so much more difficult if you want to try out the program in single or quad precision. If R8 is used consistently almost all of the work of conversion can be done simply by changing the value of the named constant R8 and recompiling. procedures that aren't generic may have to have their names changed as well to get the right specific version.

Fourth bullet point: R8 should be C_DOUBLE; both should be equal to 8. Print them out to make sure. Or you could even have defined R8 as C_DOUBLE in the first place if you wanted to.

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

First bullet point: var1=1.0E-20_R8 . You have to change the D to an E as well as appending the _R8. The decision could have been made that the kind overrode the D exponent, but it wasn't.

Second bullet point: the kind=R8 is redundant only if at least one of var1 or var2 is of type COMPLEX and at least one has kind=R8. If var1-var2 is not of type COMPLEX then REAL(var1-var2) will be default REAL: kind(REAL((var1-var2)) = kind(1.0) which probably is not what you want. If var1-var2 is of type COMPLEX, then REAL(var1-var2) just takes the real part of var1-var2, leaving its KIND intact: kind(REAL(var1-var2)) = kind(var1-var2).

Third bullet point: this is the most counterintuitive case I think. You pretty much always need the kind=R8 here, otherwise the result will be default COMPLEX: kind(cmplx(0.0_R8,0.0_R8,kind=R8)) = R8, but kind(cmplx(0.0_R8,0.0_R8)) = kind(1.0). Don't use dcmplx because it makes it so much more difficult if you want to try out the program in single or quad precision. If R8 is used consistently almost all of the work of conversion can be done simply by changing the value of the named constant R8 and recompiling. procedures that aren't generic may have to have their names changed as well to get the right specific version.

Fourth bullet point: R8 should be C_DOUBLE; both should be equal to 8. Print them out to make sure. Or you could even have defined R8 as C_DOUBLE in the first place if you wanted to.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Because you (or the compile you choose) may change the value of the r8 named constant (and you don't really care what the value is - just what that value means), avoid calling the constant r8.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Repeat Offender Thank you for your reply. I've tested the first bullet point (var=1.0E-20_R8) and it worked. I also took the suggestion of IanH and I'm now using RE8 instead of R8.

I do have another problem regarding precision when initializing variables using namelist. Since I cannot use _RE8 in the namelist (because it doesn't know what RE8 is), is there any way to force RE8 precision of the variables set in the namelist?

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

The precision on a Fortran formatted READ, which is what NAMELIST input is, is set by the variable being read into, not the string of characters being read from. Even the distinction between an 'E' exponent and a 'D' exponent is ignored. So the thing to do is to make the variables being read into have kind=RE8, and the Fortran processor should handle matters automatically from there.

BTW, I think what @IanH was driving at was that you might at some point want to switch over to single- or quad-precision variables throughout in which case a named constant containing the character '8' in it might be confusing. A convention I have seen a lot of is to have a relatively neutral name like 'wp'. Then the reader knows that he has to go back to the module where it was defined to know its meaning and won't assume that it means real variables that take up 8 bytes of storage. If you intend to stick with 8 byte real variables and at most use single- or quad-precision for testing purposes, 'R8' or 'RE8' would be fine I think.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

When I had problems finding properties including precision of various real kinds with various compilers I wrote programs called kinds.f90 (f95-compliant) and kinds03.f90 (f2003-compliant) to give all the details I could think of, for both real and integer kinds the compiler allows. You can find the programs via http://homepages.ecs.vuw.ac.nz/~harper/fortranstuff.shtml

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page