[Vm-dev] f2c/g77 and problem with FFI 64bits returning float

[Eliot Miranda]

Hi Nicolas,

    isn't the solution to simply misdeckare the fact notion as returning a double?  Provided the FFI mechanism supports all variants in use then it can collect values from all functions.  That the necessary declaration doesn't match the function' declaration when compiled with a non-conforming toolchain is unfortunate but can be lived with.

The issue is how to communicate the issue to programmers so they can diagnose and apply the work-around, which requires good documentation.

_,,,^..^,,,_ (phone)

> On Jan 26, 2017, at 4:57 PM, Nicolas Cellier <nicolas.cellier.aka.nice at gmail.com> wrote:
> 
> Hi,
> I have strange result with 64bits FFI function returning single precision float.
> Here is an example:
> 
> (LapackSGEMatrix rows: #((2.3))) absMax.
> 
> This matrix has a single element, 2.3 rounded to single precision float
> (2.299999952316284 when printed as a double precision)
> 
> absMax is supposed to take the maximum of absolute values in the matrix.
> It does so thru Lapack function slange:
> "
> *  Purpose
> *  =======
> *  SLANGE  returns the value of the one norm,  or the Frobenius norm, or
> *  the  infinity norm,  or the  element of  largest absolute value  of a
> *  real matrix A.
> "
>     <cdecl: float 'slange_'( char * long * long * float * long * float * long )>
> 
> Unfortunately above snippet returns 3.6893488147419103e19
> 
> It correctly calls this:
>             floatRet = dispatchFunctionPointerwithwithwithwithwithwith(((float (*)(sqIntptr_t, sqIntptr_t, sqIntptr_t, sqIntptr_t, sqIntptr_t, sqIntptr_t)) procAddr), ((calloutState->integerRegisters))[0], ((calloutState->integerRegisters))[1], ((calloutState->integerRegisters))[2], ((calloutState->integerRegisters))[3], ((calloutState->integerRegisters))[4], ((calloutState->integerRegisters))[5]);
> 
> which translates into something like:
> 
>     0x10833c537 <+2615>: movq   -0x28(%rbp), %rax
>     0x10833c53b <+2619>: movq   -0xe8(%rbp), %rcx
>     0x10833c542 <+2626>: movq   0xd8(%rcx), %rdi
>     0x10833c549 <+2633>: movq   -0xe8(%rbp), %rcx
>     0x10833c550 <+2640>: movq   0xe0(%rcx), %rsi
>     0x10833c557 <+2647>: movq   -0xe8(%rbp), %rcx
>     0x10833c55e <+2654>: movq   0xe8(%rcx), %rdx
>     0x10833c565 <+2661>: movq   -0xe8(%rbp), %rcx
>     0x10833c56c <+2668>: movq   0xf0(%rcx), %rcx
>     0x10833c573 <+2675>: movq   -0xe8(%rbp), %r8
>     0x10833c57a <+2682>: movq   0xf8(%r8), %r8
>     0x10833c581 <+2689>: movq   -0xe8(%rbp), %r9
>     0x10833c588 <+2696>: movq   0x100(%r9), %r9
> ->  0x10833c58f <+2703>: callq  *%rax
>     0x10833c591 <+2705>: cvtss2sd %xmm0, %xmm0
>     0x10833c595 <+2709>: movsd  %xmm0, -0x150(%rbp)
> 
> If I print $xmm0 just after the callq, then
> (lldb) nexti
> (lldb) print $xmm0
> (unsigned char __attribute__((ext_vector_type(16)))) $212 = (0x00, 0x00, 0x00, 0x60, 0x66, 0x66, 0x02, 0x40, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00)
> 
> and just after the connversion to double precision:
> (lldb) nexti
> (lldb) print $xmm0
> (unsigned char __attribute__((ext_vector_type(16)))) $213 = (0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00)
> 
> Let's see:
> 
> tmp := #[16r00   16r00   16r00   16r60   16r66   16r66   16r02   16r40 ].
> {tmp doubleAt: 1.
> tmp floatAt: 1}.
>  #(2.299999952316284 3.6893488147419103e19)
> 
> Bingo! that means that the value returned in xmm0 was already in double precision.
> When we convert it back to single precision (it's like interpreting the 4 LSB of the double as a single precision), then we get the incorrect value...
> 
> So why was slange result promoted to double?
> I can reproduce on macosx with pre-installed veclib, and in win64 compiling LAPACK 3.3.1 from sources (translated by f2c) with MSVC10.
> 
> Ah, Ah, f2c! Dont you promote float return values to double? YES
> But why this does not happen with the 32bits VM ???
> That's what drove me off the solution for a while...
> It's the IA32 ABI... return value is stored in ST0 (allways promoted to double).
> So converting it to a double again like we do is a no-op and just works in 32bits.
> 
> That's going to be a problem for FORTRAN functions on 64bits.
> IF compiled thru g77 or f2c conventions, then float results are promoted to double!
> IF compiled thru gfortran, then float result just remain float results.
> It means a major source of incompatibility: how to guess how this binary was compiled? (for example vecLib...)
> 
> And how to adapt my FFI source code?
> Last thing, f2c might also be non standard when returning a complex value
> Big ball of mud...
>