Julia and Fortran: passing floating point numbers

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Julia and Fortran: passing floating point numbers

Alexey Cherkaev
I have encountered a rather strange error as I'm trying to create a wrapper for Fortran's TWPBVP (ODE BVP) solver.

Full code is available at https://github.com/mobius-eng/TWPBVP.jl

After downloading (to julia's package directory), run

Pkg.build("TWPBVP")

and

Pkg.test("TWPBVP")

The issue (right now) is the parameter "aright" or "right": in Julia code it is set to "float(pi)", but when it is passed to Fortran's code, it becomes "0.0"

To cut long story short, I have a problem of passing floating-point numbers from Julia to Fortran subroutine. Documentation suggests, that each parameter to Fortran's program needs "Ref{T}" wrapper as they are passed by reference (I've seen that older code also uses Ptr{T} with "&" at variable names). It works fine for integers (I even can see the difference if "FInt" alias is changed between "Int32" and "Int64" and if "-fdefault-integer-8" in compiling option is present or omitted). But I don't see a single floating point value passed correctly (including floating-point arrays). I've also tried, as an extreme measure, to pass a floating point parameter as a vector (of size 1) and have a corresponding type in ccall to be "Ptr{Float64}"

SysInfo: Ubuntu 64-bit, Julia 0.5

Reply | Threaded
Open this post in threaded view
|

Re: Julia and Fortran: passing floating point numbers

Angel de Vicente
Hi,

Alexey Cherkaev <[hidden email]> writes:

> To cut long story short, I have a problem of passing floating-point numbers from
> Julia to Fortran subroutine. Documentation suggests, that each parameter to
> Fortran's program needs "Ref{T}" wrapper as they are passed by reference (I've
> seen that older code also uses Ptr{T} with "&" at variable names). It works fine
> for integers (I even can see the difference if "FInt" alias is changed between
> "Int32" and "Int64" and if "-fdefault-integer-8" in compiling option is present
> or omitted). But I don't see a single floating point value passed correctly
> (including floating-point arrays). I've also tried, as an extreme measure, to
> pass a floating point parameter as a vector (of size 1) and have a corresponding
> type in ccall to be "Ptr{Float64}"

Perhaps in Fortran you have 32bits floats? I just tried with this, and I
have no problem passing floats to Fortran:

,----
| testdp(ar) = ccall((:__tests_MOD_testdp, "./testmodgfort"),Float64, (Ptr{Float64},), &ar)
| println("Result of function call to Fortran ", testdp(5.0))
`----

,----
| DOUBLE PRECISION FUNCTION testdp (a)
|   DOUBLE PRECISION :: a
|   testdp = 5.0*a
| END FUNCTION testdp
`----

,----
| julia> include("julia-fortran.jl")
| Result of function call to Fortran 25.0
`----

--
Ángel de Vicente
http://www.iac.es/galeria/angelv/         
Reply | Threaded
Open this post in threaded view
|

Re: Julia and Fortran: passing floating point numbers

Alexey Cherkaev
On a smaller example I can see that it is indeed is because floats are not double floats.

It's seem to be the statement

implicit double precision (a-h,o-z)

in Fortran code that throws things upside down. I'm not so familiar with the Fortran standard and especially with gfortran extensions, but I would think if anything, this statement would force the arguments to be double float (real*8) precision. Also, I put "-fdefault-real-8" option to the compiler...

But thanks for the tip, perhaps I would just comment this statement out in Fortran.


On Wednesday, November 9, 2016 at 3:49:55 PM UTC+2, Ángel de Vicente wrote:
Hi,

Alexey Cherkaev <<a href="javascript:" target="_blank" gdf-obfuscated-mailto="hfZ1sZwJBQAJ" rel="nofollow" onmousedown="this.href=&#39;javascript:&#39;;return true;" onclick="this.href=&#39;javascript:&#39;;return true;">alexey....@...> writes:

> To cut long story short, I have a problem of passing floating-point numbers from
> Julia to Fortran subroutine. Documentation suggests, that each parameter to
> Fortran's program needs "Ref{T}" wrapper as they are passed by reference (I've
> seen that older code also uses Ptr{T} with "&" at variable names). It works fine
> for integers (I even can see the difference if "FInt" alias is changed between
> "Int32" and "Int64" and if "-fdefault-integer-8" in compiling option is present
> or omitted). But I don't see a single floating point value passed correctly
> (including floating-point arrays). I've also tried, as an extreme measure, to
> pass a floating point parameter as a vector (of size 1) and have a corresponding
> type in ccall to be "Ptr{Float64}"

Perhaps in Fortran you have 32bits floats? I just tried with this, and I
have no problem passing floats to Fortran:

,----
| testdp(ar) = ccall((:__tests_MOD_testdp, "./testmodgfort"),Float64, (Ptr{Float64},), &ar)
| println("Result of function call to Fortran ", testdp(5.0))
`----

,----
| DOUBLE PRECISION FUNCTION testdp (a)
|   DOUBLE PRECISION :: a
|   testdp = 5.0*a
| END FUNCTION testdp
`----

,----
| julia> include("julia-fortran.jl")
| Result of function call to Fortran 25.0
`----

--
Ángel de Vicente
<a href="http://www.iac.es/galeria/angelv/" target="_blank" rel="nofollow" onmousedown="this.href=&#39;http://www.google.com/url?q\x3dhttp%3A%2F%2Fwww.iac.es%2Fgaleria%2Fangelv%2F\x26sa\x3dD\x26sntz\x3d1\x26usg\x3dAFQjCNHWgWJwB7xZkScAhM-XaTehlCOVww&#39;;return true;" onclick="this.href=&#39;http://www.google.com/url?q\x3dhttp%3A%2F%2Fwww.iac.es%2Fgaleria%2Fangelv%2F\x26sa\x3dD\x26sntz\x3d1\x26usg\x3dAFQjCNHWgWJwB7xZkScAhM-XaTehlCOVww&#39;;return true;">http://www.iac.es/galeria/angelv/          
Reply | Threaded
Open this post in threaded view
|

Re: Julia and Fortran: passing floating point numbers

Alexey Cherkaev
Found the problem: "double precision" on my machine (presumably, any 64-bit machine) with "-fdefault-real-8" flag will be "real*16" (something like non-existent "Float128" in Julia).

On Wednesday, November 9, 2016 at 4:05:59 PM UTC+2, Alexey Cherkaev wrote:
On a smaller example I can see that it is indeed is because floats are not double floats.

It's seem to be the statement

implicit double precision (a-h,o-z)

in Fortran code that throws things upside down. I'm not so familiar with the Fortran standard and especially with gfortran extensions, but I would think if anything, this statement would force the arguments to be double float (real*8) precision. Also, I put "-fdefault-real-8" option to the compiler...

But thanks for the tip, perhaps I would just comment this statement out in Fortran.


On Wednesday, November 9, 2016 at 3:49:55 PM UTC+2, Ángel de Vicente wrote:
Hi,

Alexey Cherkaev <[hidden email]> writes:

> To cut long story short, I have a problem of passing floating-point numbers from
> Julia to Fortran subroutine. Documentation suggests, that each parameter to
> Fortran's program needs "Ref{T}" wrapper as they are passed by reference (I've
> seen that older code also uses Ptr{T} with "&" at variable names). It works fine
> for integers (I even can see the difference if "FInt" alias is changed between
> "Int32" and "Int64" and if "-fdefault-integer-8" in compiling option is present
> or omitted). But I don't see a single floating point value passed correctly
> (including floating-point arrays). I've also tried, as an extreme measure, to
> pass a floating point parameter as a vector (of size 1) and have a corresponding
> type in ccall to be "Ptr{Float64}"

Perhaps in Fortran you have 32bits floats? I just tried with this, and I
have no problem passing floats to Fortran:

,----
| testdp(ar) = ccall((:__tests_MOD_testdp, "./testmodgfort"),Float64, (Ptr{Float64},), &ar)
| println("Result of function call to Fortran ", testdp(5.0))
`----

,----
| DOUBLE PRECISION FUNCTION testdp (a)
|   DOUBLE PRECISION :: a
|   testdp = 5.0*a
| END FUNCTION testdp
`----

,----
| julia> include("julia-fortran.jl")
| Result of function call to Fortran 25.0
`----

--
Ángel de Vicente
<a href="http://www.iac.es/galeria/angelv/" rel="nofollow" target="_blank" onmousedown="this.href=&#39;http://www.google.com/url?q\x3dhttp%3A%2F%2Fwww.iac.es%2Fgaleria%2Fangelv%2F\x26sa\x3dD\x26sntz\x3d1\x26usg\x3dAFQjCNHWgWJwB7xZkScAhM-XaTehlCOVww&#39;;return true;" onclick="this.href=&#39;http://www.google.com/url?q\x3dhttp%3A%2F%2Fwww.iac.es%2Fgaleria%2Fangelv%2F\x26sa\x3dD\x26sntz\x3d1\x26usg\x3dAFQjCNHWgWJwB7xZkScAhM-XaTehlCOVww&#39;;return true;">http://www.iac.es/galeria/angelv/          
Reply | Threaded
Open this post in threaded view
|

Re: Julia and Fortran: passing floating point numbers

Angel de Vicente
Alexey Cherkaev <[hidden email]> writes:

> Found the problem: "double precision" on my machine (presumably, any 64-bit
> machine) with "-fdefault-real-8" flag will be "real*16" (something like
> non-existent "Float128" in Julia).

In Fortran, if you want your code to be portable and avoid these type of
issues it is much better to use selected_real_kind (see
https://gcc.gnu.org/onlinedocs/gfortran/SELECTED_005fREAL_005fKIND.html)

Cheers,
--
Ángel de Vicente
http://www.iac.es/galeria/angelv/