REAL (KIND= )

[reidar]

In a calculation procedure I usually apply REAL*8  in math operation. However sometimes accuracy/round off errors affect the results. Then I discovered that with IVF REAL*16 is accepted.  This brings accuracy to a new level, but the downside  is of course that computation time is tripled, at least!

What I wander is, if it is possible to set the KIND=8 or 16 at the time of computation, as a part of the input data?   

I made  test with REAL (KIND=NOTSET) retset

error #6683: A kind type parameter must be a compile-time constant. [NOTSET]

I think this indicates that the KIND parameter must be set to 4, 8 or 16 before compilation. Or is there a way around?

Regards

[Arjen_Markus]

No, this is a compile-time constant, not something you can change on the fly.

However, it is possible to provide alternative implementations from the same code, with merely a different kind, and select the one you want via a SELECT CASE or the like. It is a trifle involved though.

[reidar]

Thank you for your response Markus Arjen. However, I du not understand how a SELECT CASE can solve it without having two different versions of most variables and routines. One set with 8 bytes and one with 16 bytes variables. 

[Arjen_Markus]

Something along these lines:

real(8) ::x , y
select case (prec)
     case(8)
           x = f(y)
     case(16)
           x = f( real(y,kind=16))
end select

where f is an interface to the two underlying implementations

[JVanB]

To create two different versions of subprograms you need some form of template capability in the language. Fortran doesn't have an explicit capability, but it can be faked to a certain extent via modules. Subprograms that require the two forms are modified so that the KIND is given by a named constant that will be acquired via host association and then they are place in INCLUDE files. Modules are written which just set the named constant to the appropriate value and include the files from the previous sentence. Then you just need a root subprogram which your main program can invoke that starts working with the right KINDs. Let's see... do I have time to compose a simple example?

Suppose your program looked something like this:

! hello1.f90
program hello1
   use ISO_FORTRAN_ENV, only: dp => REAL64
   implicit none
   real(dp) pi
   pi = 4*atan(1.0_dp)
   write(*,*) 'Hello, world: pi = ',pi
end program hello1

First we put the good stuff in an INCLUDE file
 

! hello2.i90
subroutine hello
   real(wp) pi
   pi = 4*atan(1.0_wp)
   write(*,*) 'Hello, world: pi = ',pi
end subroutine hello

Note that we have changed the name of the kind type parameter (not really necessary) and used the name consistently throughout. This is our template version of the program.

Then we have to tell the compiler to make a version for each KIND we want and then create a module that renames the root subprogram and then our program can decide which one to invoke.

! hello3.f90
module M8
   use ISO_FORTRAN_ENV, only: wp => REAL64
   implicit none
   contains
include 'hello2.i90'
end module M8

module M16
   use ISO_FORTRAN_ENV, only: wp => REAL128
   implicit none
   contains
include 'hello2.i90'
end module M16

module M
   use M8, only: hello8 => hello
   use M16, only: hello16 => hello
   implicit none
end module M

program P
   use M
   implicit none
   integer choice
   write(*,'(a)',advance='no') 'Enter 2 for quad precision:> '
   read(*,*) choice
   if(choice == 2) then
      call hello16
   else
      call hello8
   end if
end program P

Of course there are many possible variations of this technique but hopefully this can get you started in the right direction.

[reidar]

Thank you for response, Repeat Offender. I  think you want to give a hint about how to set default real variables to either 4, 8 or 16 bytes?  As I am not familiar to the module ISO_Fortran, I just added a use statement to a random file as shown below. 

    USE ISO_FORTRAN_ENV  , only: wp -> REAL64

I received the following response:    
    C:\CALC\CaFeMS\Main_elast6b.f90(68): error #5082: Syntax error, found '-' when expecting one of: ( , <END-OF-STATEMENT> ; (/ =>
compilation aborted for C:\CALC\CaFeMS\Main_elast6b.f90 (code 1)
compilation aborted for C:\CALC\CaFeMS\Main_elast6b.f90 (code 1)

Hope you can give me some more clues.

[Arjen_Markus]

It is a simple syntax error:

only: wp -> REAL64

shoudl be:

only: wp => REAL64

[reidar]

Repeat Offender: I understand you try to show how to extract a global parameter wp (module M8 and M16) that can be contained by another module (?) that then again is assigned in the actual subroutines. Right ?

The assignment of this global parameter takes place in hello8/hello16 subroutines.  But I do not understand the link between the (undefined) hello, hello8 and hello16 routines.. or are these "routines" just symbol variables that invoke different parts in the M module ?

[JVanB]

OK reidar, let's try reasoning this out step by step.  First, after the INCLUDEs are processed by the compiler, it's going to be working on something like this:

! hello3.f90
module M8
   use ISO_FORTRAN_ENV, only: wp => REAL64
   implicit none
   contains
!include 'hello2.i90'
! hello2.i90
subroutine hello
   real(wp) pi
   pi = 4*atan(1.0_wp)
   write(*,*) 'Hello, world: pi = ',pi
end subroutine hello
end module M8

module M16
   use ISO_FORTRAN_ENV, only: wp => REAL128
   implicit none
   contains
!include 'hello2.i90'
! hello2.i90
subroutine hello
   real(wp) pi
   pi = 4*atan(1.0_wp)
   write(*,*) 'Hello, world: pi = ',pi
end subroutine hello
end module M16

module M
   use M8, only: hello8 => hello
   use M16, only: hello16 => hello
   implicit none
end module M

program P
   use M
   implicit none
   integer choice
   write(*,'(a)',advance='no') 'Enter 2 for quad precision:> '
   read(*,*) choice
   if(choice == 2) then
      call hello16
   else
      call hello8
   end if
end program P

Notice that the INCLUDEs have been replaced by copies of the file that was named in the INCLUDE statements. Thus both module M8 and M16 will posses a module procedure called SUBROUTINE HELLO. There isn't really a problem of ambiguity, however, because the compiler will mangle their names to something like M8_mp_HELLO and M16_mp_HELLO. They must have the same name as seen by Fortran within their modules, however, because they were sourced from the same INCLUDE file. If they are going to be visible and used in the same scope (here PROGRAM P) one or both must be renamed. That is the job of MODULE M: it provides names visible to any program units that USE it, HELLO8 for the double precision version that was compiled in MODULE M8 and HELLO16 for the quad precision version from MODULE M16. Both of these modules had a procedure that had the same name (WP) for the REAL KIND we wanted for that procedure, but the named constant picked up its value via host association from its module and was REAL64 [=8] from the specification part of M8 and REAL128 [=16] from the specification part of M16.

So now if we USE M, we have two different names accessible to us: HELLO8 which really means the SUBROUTINE HELLO that was compiled in MODULE M8 where wp had the value 8, and HELLO16 which refers to the SUBROUTINE HELLO that was compiled in MODULE M16 where wp had the value 16. You might be thinking that something more profound is going on here, like C++ templates, but if you think about it carefully you can see how we end up with two different subroutines, HELLO8 and HELLO16, visible in PROGRAM P which were created from the same source code, but associated with different hosts.

[jimdempseyatthecove]

For function and subroutines, you can also declare a generic interface, say

interface f
  function f4(x) result(y)
    real(4), intent(in) :: x
    real(4) :: y
  end function f4
  function f8(x) result(y)
    real(8), intent(in) :: x
    real(8) :: y
  end function f8
  function f16(x) result(y)
    real(16), intent(in) :: x
    real(16) :: y
  end function f16
end interface

You will need to write the appropriate routines (possibly using the INCLUDE to insert common code).

Generic interface is an exposed means to perform function overloading such as in C/C++/C#.

Jim Dempsey

[reidar]

Thanks for comments Jimdempsey, I will study.

But further to advise from Repeat Offender: I copied and compiled/run your first #5 post. Post #9 is similar but now the hello routine is fully written in the module. I think I understand how it works. But how do I transfer the parameter "wp" to all the other subroutines to declare variables as REAL*8 or REAL *16, so that I can declare the actual  variables with  REAL (KIND=wp) ?  Suppose the "hello" subroutine is number- crunching program that calls other functions and subroutines in which the selected variables are to be declared as 8 or 16 bytes variables.   Is that possible or do I have to arrange assignment of the M8 alternatively M16 module in each subroutine?

[jimdempseyatthecove]

Repeat Offender's suggestion showed you how to create a facsimile of  a template for use as an INCLUDE file used in a CONTAINS section of a MODULE qualified by a USE ISO_FORTRAN_ENV, only: wp=>...yourDesiredPrecisionHere...

If you want to globally specify a precision, then one of the many ways to do this would be to

a) choose what you want for an identifier, I will choose "wp", you can choose what you want
b) write a module file, perhaps you can call it defaults.f90, and in there have in the data (pre-contains) section
      USE ISO_FORTRAN_ENV, only: wp=>...yourDesiredPrecisionHere...
c) then add "USE defaults" to the subroutines and functions.
d) I assume you understand by now that the module file will have to be compiled and be a dependency of your program.

Note, you likely have a common module that you already USE. You can insert the "USE defaults" inside that if you wish, or simply insert  USE ISO_FORTRAN_ENV, only: wp=>...yourDesiredPrecisionHere...

Jim Dempsey

[JVanB]

I tried to make an example that does what you want, but encountered problems possibly due to compiler bugs. Here is what I did: start out with a program that does a little calculation:

subroutine cubic(a,b,c,d,x,n)
   implicit none
   real a,b,c,d
   real x(*)
   integer n
! Solves a*x**3+b*x**2+c*x+d = 0
! Assumes a /= 0
! x(1:n) are the n real solutions, sorted ascending
   real p, q, offset
   real r
   real, parameter :: pi = 4*atan(1.0)

   offset = b/(3*a)
   p = c/a-b**2/(3*a**2)
   q = 2*b**3/(27*a**3)-b*c/(3*a**2)+d/a
   if(p > 0) then
      r = sqrt(4*p/3)
      n = 1
      x(1:n) = r*sinh(asinh(-4*q/r**3)/3)-offset
   else if(p == 0) then
      n = 1
      x(1:n) = -q**(1.0/3)-offset
   else ! p < 0
      r = sign(sqrt(-4*p/3),-q)
      if(q == 0) then
         n = 3
         x(1:n) = [-sqrt(-p),0.0,sqrt(-p)]-offset
      else if(abs(4*q) < abs(r**3)) then
         n = 3
         x(1:n) = r*cos(acos(-4*q/r**3)/3+[2*pi/3,-2*pi/3,0.0])-offset
         if(r < 0) x(1:n) = x(n:1:-1)
      else if(abs(4*q) == abs(r**3)) then
         n = 2
         x(1:n) = r*[-0.5,1.0]-offset
         if(r < 0) x(1:n) = x(n:1:-1)
      else ! abs(4*q) > abs(r**3)
         n = 1
         x(1:n) = r*cosh(acosh(-4*q/r**3)/3)-offset
      end if 
   end if
end subroutine cubic

module vdw
   implicit none
   type params
      character(:), allocatable :: Name
      real a ! kPa**L**2/mol**2
      real b ! L/mol
   end type params
   contains
      function v(P,T,ab)
         real P, T
         type(params) ab
         real, parameter :: R = 8.3144598 ! kPa*L/(K*mol)
         real v
         real x(3)
         integer n
         call cubic(P,-(ab%b*P+R*T),ab%a,-ab%a*ab%b,x,n)
         v = x(n)
      end function v
end module vdw

program vdw1
   use vdw
   implicit none
   real P
   real T
   type(params) ab
   real molar_volume

   write(*,'(a)',advance='no') 'Enter P(kPA) :> '
   read(*,*) P
   write(*,'(a)',advance='no') 'Enter T(K) :> '
   read(*,*) T
   ab = params('Xe',425.0,0.05105)
   molar_volume = v(P,T,ab)
   write(*,*) ab%name,molar_volume,P,8.3144598*T/ &
      (molar_volume-ab%b)-ab%a/molar_volume**2
end program vdw1

I don't know if I'm doing the right thing here: when van der Waals' equation spits out multiple roots for the molar volume are you supposed to take the maximum root? Doesn't really matter as this is just an example. So if you compile the above and enter a reasonable pressure and temperature at the prompts, it seems to yield at least a consistent molar volume for xenon. Again, the first step is to convert all REAL KINDs to a named constant. If all REAL literals have a decimal point it's much easier to search them out and append the _wp in each case.

module default
   use ISO_FORTRAN_ENV, only: wp => REAL64
end module default

subroutine cubic(a,b,c,d,x,n)
   use default
   implicit none
   real(wp) a,b,c,d
   real(wp) x(*)
   integer n
! Solves a*x**3+b*x**2+c*x+d = 0
! Assumes a /= 0
! x(1:n) are the n real solutions, sorted ascending
   real(wp) p, q, offset
   real(wp) r
   real(wp), parameter :: pi = 4*atan(1.0_wp)

   offset = b/(3*a)
   p = c/a-b**2/(3*a**2)
   q = 2*b**3/(27*a**3)-b*c/(3*a**2)+d/a
   if(p > 0) then
      r = sqrt(4*p/3)
      n = 1
      x(1:n) = r*sinh(asinh(-4*q/r**3)/3)-offset
   else if(p == 0) then
      n = 1
      x(1:n) = -q**(1.0_wp/3)-offset
   else ! p < 0
      r = sign(sqrt(-4*p/3),-q)
      if(q == 0) then
         n = 3
         x(1:n) = [-sqrt(-p),0.0_wp,sqrt(-p)]-offset
      else if(abs(4*q) < abs(r**3)) then
         n = 3
         x(1:n) = r*cos(acos(-4*q/r**3)/3+[2*pi/3,-2*pi/3,0.0_wp])-offset
         if(r < 0) x(1:n) = x(n:1:-1)
      else if(abs(4*q) == abs(r**3)) then
         n = 2
         x(1:n) = r*[-0.5_wp,1.0_wp]-offset
         if(r < 0) x(1:n) = x(n:1:-1)
      else ! abs(4*q) > abs(r**3)
         n = 1
         x(1:n) = r*cosh(acosh(-4*q/r**3)/3)-offset
      end if 
   end if
end subroutine cubic

module vdw
   use default
   implicit none
   type params
      character(:), allocatable :: Name
      real(wp) a ! kPa**L**2/mol**2
      real(wp) b ! L/mol
   end type params
   contains
      function v(P,T,ab)
         real(wp) P, T
         type(params) ab
         real(wp), parameter :: R = 8.3144598_wp ! kPa*L/(K*mol)
         real(wp) v
         real(wp) x(3)
         integer n
         call cubic(P,-(ab%b*P+R*T),ab%a,-ab%a*ab%b,x,n)
         v = x(n)
      end function v
end module vdw

program vdw2
   use default
   use vdw
   implicit none
   real(wp) P
   real(wp) T
   type(params) ab
   real(wp) molar_volume

   write(*,'(a)',advance='no') 'Enter P(kPA) :> '
   read(*,*) P
   write(*,'(a)',advance='no') 'Enter T(K) :> '
   read(*,*) T
   ab = params('Xe',425.0_wp,0.05105_wp)
   molar_volume = v(P,T,ab)
   write(*,*) ab%name,molar_volume,P,8.3144598_wp*T/ &
      (molar_volume-ab%b)-ab%a/molar_volume**2
end program vdw2

The problem with this approach is that you have to recompile everything to convert to another KIND, although the only change that need be made is in the second line where the wp is determined and then propagated throughout the program. So I wanted to make an example that doesn't duplicate code but doesn't require all the calculational stuff to be in an INCLUDE file. The idea is that I can write out a file default8.f90 that has that module default in it that sets wp => REAL64. After compiling this we can compile calc.f90 to calc8.dll. Since it USEes the module default that we just compiled, its module vdw will have its wp = REAL64 and its type(params) will have REAL(REAL64) components and so on. Then we compile types8.f90 which gives us a module types8 which is really an alias for the module vdw with wp =" REAL64 inside.

OK, then we repeat the procedure, compiling default16.f90 so now we will have a module default with wp => REAL128. Now when we compile calc.f90 to calc16.dll we get a module vdw with wp=> REAL128 inside. Then when we compile types16.f90 we get a module types16 that is an alias for the most recent modulw vdw we compiled, with wp => REAL128 internally. Thus our only type-specific code to this point was in the 4 files default8.f90, types8.f90, default16.f90, and types16.f90, 3 lines each. This is good from the standpoint of code maintenance.

Now our main program is written to use module types8, and we use LoadLibrary and GetProcAddress to get a pointer to the entry point of calc8.dll. When the user responds to the prompt saying that he wants to do the quad precision calculation, the BLOCK construct is entered that gets type info from module types16 and now gets a pointer to the entry point of calc16.dll so now it can just do the computational stuff that was compiled in quad precision. A little more boilerplate in the main program and we need those 4 extra little files, but the part about cahnging the big computational part to INCLUDE files is gone because there are no INCLUDE file any more.

Unfortunately both gfortran and ifort die on this example. When gfortran is invoked via the batch file gf_vdw.bat, the first and only error I get is

vdw3.f90:63: confused by earlier errors, bailing out

and when ifort is invoked via make_vdw3.bat, the error parade starts with

vdw3.f90(41): error #6405: The same named entity from different modules and/or p
rogram units cannot be referenced.   [WP]
   real(wp) P
--------^

Here is a *.zip file with the failing example:

[jimdempseyatthecove]

Try commenting out line 70 in program (use default is used in vdw)

This shouldn't have been an error unless the date and time on the .mod file used to build the library differed from the date and time used to build the program.

Jim Dempsey

[JVanB]

I'm not sure what you mean about line 70. Line 70 of vdw3.f90 is the first half of a WRITE statement program vdw2 is exactly 70 lines long so this would be the end line.

I think what's happening is that ifort doesn't build a *.mod file containing full info about all symbols but rather it only provides references for symbols it acquires via USE association. And you can't have two different *.mod files with the same name in the visible INCLUDE directories. I tried changing make_vdw3.bat to

ifort /nologo /c default8.f90
IF EXIST sd8\NUL GOTO HAVEsd8
mkdir sd8
:HAVEsd8
move default.mod sd8
ifort /nologo /dll /Isd8 calc.f90 /exe:calc8
move vdw.mod sd8
ifort /nologo /c /Isd8 types8.f90
move types8.mod sd8
ifort /nologo /c default16.f90
IF EXIST sd16\NUL GOTO HAVEsd16
mkdir sd16
:HAVEsd16
move default.mod sd16
ifort /nologo /dll /Isd16 calc.f90 /exe:calc16
move vdw.mod sd16
ifort /nologo /c /Isd16 types16.f90
move types16.mod sd16
ifort /nologo /Isd8 /Isd16 vdw3.f90

And now ifort dies at

vdw3.f90(77): error #6405: The same named entity from different modules and/or p
rogram units cannot be referenced.   [BLOCK_WP]
         real(BLOCK_wp) BLOCK_P
--------------^

Because the first default.mod file in its search path is sd8\default.mod so it was OK when the wp from that module file was required but now when it needs the wp from sd16\default.mod it finds the wrong one because the right one is later in the include search path even though the module that is in the chain leading to it is in sd16. So I can't see a workaround for this issue.

[jimdempseyatthecove]

>>I'm not sure what you mean about line 70. Line 70 of vdw3.f90

Perhaps I should have made this a bit more clear:

In your post #13, you have two code examples, in the second code example, you have several program units, listed in a single paste. In this second example, on line 70 (which is not necessarily the 70'th line of the program, but is likely unknown as to which line it is), you have what may be a potential incorrect USE (if the module used differs in date/time from the same named module used in the used module named vdw).

Remove/comment the "use default" from "program vdw2"

Jim Dempsey

[JVanB]

Oh, that example. No, that's not an incorrect USE because in that case there was only one instance of module default, right there at the top of the file vdw2.f90. That example compiles and runs OK, but requires recompilation to change the REAL KIND. The one the doesn't work is in the vdw3.zip attachment to the post. That one generates two versions of default.mod and there is no way of telling ifort which one to use if both are necessary for a single invocation of the compiler, either directly or chained through USE in other modules. One could write an independent module for use with the final program, but this would be like what you have to do in C where you maintain separate *.c and *.h files, but with the disadvantage that you couldn't confirm the consistency of the *.h files with the *.c files because of the generic nature of the *.c (actually *.f90) files in this case. An extra maintenance headache on top of the already somewhat more complicated build sequence.

Perhaps you could have tried compiling my example so that you better understood what you were commenting on.

[jimdempseyatthecove]

The problem may be with the "association" wp => ...

Instead of using the "association", for each compilation unit (type) use an INTEGER, PARAMETER :: wp=...

In this manner, wp should be static to the compilation unit (together with the modules it uses in the compilation).

Jim Dempsey

[Steve_Lionel]

Repeat Offender wrote:

That one generates two versions of default.mod and there is no way of telling ifort which one to use if both are necessary for a single invocation of the compiler, either directly or chained through USE in other modules. 

Maybe you know this, but the language disallows having more than one program unit (a module is a program unit) of the same name in a program. (F2008 16.2). Note that an intrinsic module isn't a program unit, so it's ok to have your own module with the same name as an intrinsic module. I would argue that compiling the same source more than once in a build constitutes a duplication as far as this is concerned.

[JVanB]

@jimdempseyatthecove, the problem is that the object of the exercise is to reuse the same file to compile a different function, so the rules of the game dictate that every statement be the same.

@Steve Lionel (Ret.), of course you know that only makes my determination greater! I did get something to work in vdw4.zip, attached below. To build vdw4.exe, just run make_vdw4.bat. Some notes:

I had a bug where I computed the cube root as

x(1:n) = -q**(1.0_wp/3)-offset

When I should have written

x(1:n) = sign(abs(q)**(1.0_wp/3),-q)-offset

but the surprising thing was that ifort interpreted my original code as a request to compute the cube root, so my cubic equation solver didn't crash even when p == 0 and q < 0 :)

Also I forgot to NUL-terminate those strings passed to LoadLibrary and GetProcAddress.

I had the wrong variable P instead of BLOCK_P in my last write statement.

I fixed these issues and restructured as threatened in my last post so that new interfaces are written ab initio with the help of an INCLUDE file.

calc.f90 is still compiled twice in this build.

One might be able to claim that it's OK to have two program units with the same name in two different *.DLLs even if the final program is going to dynamically link to both of them. Has the latest standard gotten around to allowing non-BIND(C) UDTs as dummy arguments to BIND(C) procedures? That would make it easier to rename them without !DEC$ ATTRIBUTES ALIAS. The usage with these UDTs really is nonstandard, because the compiler could reorder their components. You can't have a SEQUENCE type with an allocatable component, can you?

[jimdempseyatthecove]

module default32
   use ISO_FORTRAN_ENV, only: REAL32
   integer, parameter :: wp = REAL32
end module default32
    
module default64
   use ISO_FORTRAN_ENV, only: REAL64
   integer, parameter :: wp = REAL64
end module default64
    
module default128
   use ISO_FORTRAN_ENV, only: REAL128
   integer, parameter :: wp = REAL128
end module default128

! ================== different source file ==========
module default
    use default64
end module default
    
! ================== different source file ==========
program wptest
    use default
    implicit none
    real(wp) :: x
    x = 123.4_wp
end program wptest

Try revising your code to use something like the above.

Jim Dempsey