R problems with lapack with gfortran

[Thomas König]

Hi,

I have tried to pinpoint potential problems which could lead to the
LAPACK issues that are currently seen in R.  I built the current R
trunk using

AR=gcc-ar RANLIB=gcc-ranlib ./configure --prefix=$HOME --enable-lto
--enable-BLAS-shlib=no --without-recommended-packages

and used this to find problem areas.

There are quite a few warnings that were flagged, due to mismatches
in function types.

The prototypes that R has in its header files, for example BLAS.h,
are often not compatible with gfortran function declarations.  To take
one small example, in src/main/print.c, we have

void NORET F77_NAME(xerbla)(const char *srname, int *info)

so xerbla_ is defined with two arguments.

However, gfortran passes string lengths as hidden arguments.
You can see this by compiling the small example

$ cat xer.f
       SUBROUTINE FOO
       INTEGER INFO
       CALL XERBLA ('FOO', INFO)
       END
$ gfortran -c -fdump-tree-original xer.f
$ cat xer.f.004t.original
foo ()
{
   integer(kind=4) info;

   xerbla (&"FOO"[1]{lb: 1 sz: 1}, &info, 3);
}

so here we have three arguments. This mismatch is flagged
by -Wlto-type-mismatch, which, for example, yields

print.c:1120:12: note: type 'void' should match type 'long int'
../../src/extra/blas/blas.f:357:20: warning: type of 'xerbla' does not
match original declaration [-Wlto-type-mismatch]
   357 |          CALL XERBLA( 'DGBMV ', INFO )


So, why can gcc's r268992 / r269349 matter? Before these patches,
gfortran used the variadic calling convention for calling procedures
outside the current file, and the non-variadic calling convention for
calling procedures found in the current file.

Because the procedures were all compiled as non-variadic, the caller and
the calle's signature did not match if they were not in the same
source file, which is an ABI violation.

This violation manifested itself in https://gcc.gnu.org/PR87689 ,
where the the problem resulted in crashes on a primary gcc platform,
POWER.

How can this potentially affect R?  After the fix for PR87689,
gfortran's calls to external procedures are no longer variadic.  It is
quite possible that, while this "works" most of the time, there
is a problem with a particular LAPACK routine, the call sequence
leading up to it or the procedures it calls.

How to fix this problem?  The only clear way I see is to fix this
on the R side, by adding the string lengths to the prototypes.
These are size_t (64 bit on 64-bit systems, 32 bit on 32-bit
systems).  You should then try to make --enable-lto pass
without any warnings.

Regarding LAPACK itself, the default build system for R builds
it as a shared library.  Offhand, I did not see any way to
build a *.a file instead, so I could not use LTO to check
for mismatched prototypes between R and LAPACK.

Of course, I cannot be sure that this is really the root cause
of the problem you are seeing,but it does seem to fit quite well.
I hope this analysis helps in resolving this.

Regards

        Thomas

[Tomas Kalibera-2]

Dear Thomas,

thank you for your input. I've debugged one of the packages and I
confirm that the breakage is related to passing of strings from C to
Fortran. Indeed, BLAS and LAPACK define a large number of subroutines
that take one or more explicit single-character strings as arguments.
Other than that, BLAS has only one function (xerbla), which takes a
string of unspecified length, LAPACK only has four (ilaenv,
ilaenv2stage, lsamen, xerbla). The C interfaces to BLAS/LAPACK from
Netlib depend on the historic behavior that explicit single-character
strings are interoperable, concretely CBLAS and LAPACKE provide C
interfaces/code that calls into Fortran BLAS/LAPACK without passing the
1s as lengths of the character strings (functions taking a string of
unspecified length are trivial and re-implemented in C). This has been
working fine for very many years as the Fortran code never needed to
access the length it knew was 1. R has been using the same practice,
which long predates ISO_C_BINDING/BIND(C), and I've seen online
discussions where people assumed interoperability of length 1 strings,
once mentioning also a citation from Fortran 2003 Handbook that says "A
Fortran character string with a length other than 1 is not
interoperable" (which invites interpretation that length 1 strings were
). I am not an expert to say whether the current Fortran standard
requires that interoperability and I assume that it does not given this
gfortran change.

This gfortran change breaks this interoperability: if a C function calls
a Fortran function, passing it a single-character string for a parameter
taking explicit single-character Fortran string, it may crash. I've
debugged one case with R package BDgraph, this example
"library(BDgraph); data.sim <- bdgraph.sim( n = 70, p = 5, size = 7, vis
= TRUE )" crashes due to corruption of C stack by Fortran function
DPOSV, when compiled with the new gfortran and with -O2. To see the
problem, one can just look at the disassembly of DPOSV (LAPACK), neither
the package nor R is not necessary:

SUBROUTINE DPOSV( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
CHARACTER          UPLO

In one case, DPOSV calls DPOTRS before returning. The new gfortran with
-O2 performs tail-call optimization, jumping to DPOTRS. In the annotated
disassembly snippet, at 11747f1, DPOSV tries to ensure that there is
constant 1 as string length of UPLO when tail-calling into DPOTRS, so it
writes it to stack where there already should have been 1 as length of
UPLO passed to DPOSV. But this argument has not been passed to DPOSV, so
this causes stack corruption.

  1174ce:       0f 85 62 ff ff ff       jne    117436 <dposv_+0xb6> <== jump if ERROR

          CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )

   1174d4:       48 8b 04 24             mov    (%rsp),%rax <======= rax holds LDB

   1174d8:       4c 89 7c 24 68          mov    %r15,0x68(%rsp) <=== save INFO to output param

   1174dd:       49 89 d8                mov    %rbx,%r8 <========== pass LDA as LDA

   1174e0:       4c 89 e1                mov    %r12,%rcx <========= pass A as A

   1174e3:       4c 8b 4c 24 08          mov    0x8(%rsp),%r9 <===== pass B as B

   1174e8:       4c 89 ea                mov    %r13,%rdx <========= pass NRHS as NRHS

   1174eb:       48 89 ee                mov    %rbp,%rsi <========= pass N as N

   1174ee:       4c 89 f7                mov    %r14,%rdi <========= pass UPLO as UPLO

   1174f1:       48 c7 44 24 70 01 00    movq   $0x1,0x70(%rsp) <=== pass 1 hidden arg on stack CORRUPTS C STACK

   1174f8:       00 00

   1174fa:       48 89 44 24 60          mov    %rax,0x60(%rsp) <=== pass LDB as LDB (stack)

       END

   1174ff:       48 83 c4 28             add    $0x28,%rsp <== remove 5 vars from stack (sframe)

   117503:       5b                      pop    %rbx

   117504:       5d                      pop    %rbp

   117505:       41 5c                   pop    %r12

   117507:       41 5d                   pop    %r13

   117509:       41 5e                   pop    %r14

   11750b:       41 5f                   pop    %r15 <=== restore register to level before call

          CALL DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )

   11750d:       e9 de 56 ef ff          jmpq   cbf0 <dpotrs_@plt> <=== tail call to dpotrs

Note that DPOSV never uses the length of the string (UPLO) from the
hidden argument, the compiler clearly knows that its length is 1. In
calls where the length is passed in registers, this does not cause
trouble (like LSAME) and indeed is needed as the registers have
different values

       IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN

   117448:       b9 01 00 00 00          mov    $0x1,%ecx

   11744d:       ba 01 00 00 00          mov    $0x1,%edx

   117452:       48 8d 35 bb 12 09 00    lea    0x912bb(%rip),%rsi        # 1a8714 <ipivot.4261+0xd14>

   117459:       4c 89 f7                mov    %r14,%rdi

   11745c:       e8 1f 3d ef ff          callq  b180 <lsame_@plt>

but it seems to me that the compiler could just refrain from setting the
length to be 1 on the stack at 1174f1, since it knows it should have
already been there. It would be a nice property if Fortran code that
never accesses the hidden arguments with the lengths of the strings,
because it knows what those lengths are, would also never write to those
hidden arguments on the stack when it knows what they are (should be).

Before the gfortran change, DPOSV would call to DPOTRS normally (no
tail-call optimization), so this problem would not occur (I tested with
268974). By disabling tail call optimization via
-fno-optimize-sibling-calls, the problem goes away also for other
packages my colleagues have identified as crashing with the new
gfortran. Did you know of any other optimization that could break this
interoperability of 1-length strings? It would be really nice to users
if this interoperability could be preserved, and if not by default than
at least with some option.

Traditionally, BLAS/LAPACK implementations are interchangeable at
dynamic linking time, using the Fortran interface that is however also
used from C, without passing lengths for fixed 1-character strings. R
supports this too, at least on some Linux distributions including
Debian/Ubuntu it is packaged so that it runs with the BLAS/LAPACK
implementation installed on the system. Even though this is probably not
correct wrt to the todays Fortran standard (I don't know for sure), this
is the common practice, and fixing this would not be easy - one would
have to create a new interface to be used from C, separate from the
Fortran one, and all software would have to start using that interface
from C. In the current situation when the Fortran interface is used,
confusion will arise with this gfortran change as different BLAS/LAPACK
implementations are built by different Fortran compilers and use a
different mix of Fortran/C for different computational subroutines. Note
CBLAS could not be readily used as it itself breaks with the current
gfortran change as well.

The same interoperability considerations apply to R packages, which
include native code that calls from C or from Fortran into the (same)
Fortran interface of BLAS/LAPACK. There would have to be a commonly
accepted C interface instead by the BLAS/LAPACK implementations, and all
of these packages would have to be modified to use that interface. If we
created such a C interface just inside R and asked all package
maintainers to update their packages, we would still have the problem
with substitution of external BLAS(/LAPACK) implementations at dynamic
linking time.

Indeed, it would be very hard to identify these problems by testing,
because at least now the crashes are quite rate (for the tail-call
optimization, a number of conditions have to be met to cause memory
corruption, first the tail optimization has to happen, then the number
of arguments has to be so large (on x86) that the lengths are passed on
stack and not in registers, we have to be lucky for the memory
corruption to map to a crash, etc).

So, any help we could get from you would be highly appreciated, be it
just a compile option to keep the old behavior or an assurance that we
are fine if we just disable the tail-call optimization. Appreciated by
us but I believe also many others who use or develop BLAS/LAPACK, but
may not have yet run into the problem, as they may not have been
regularly testing bleeding-edge versions of compilers or may not have
such a large code base to test as we have on CRAN.

Thanks
Tomas


On 4/24/19 11:32 PM, Thomas König wrote:

> Hi,
>
> I have tried to pinpoint potential problems which could lead to the
> LAPACK issues that are currently seen in R.  I built the current R
> trunk using
>
> AR=gcc-ar RANLIB=gcc-ranlib ./configure --prefix=$HOME --enable-lto
> --enable-BLAS-shlib=no --without-recommended-packages
>
> and used this to find problem areas.
>
> There are quite a few warnings that were flagged, due to mismatches
> in function types.
>
> The prototypes that R has in its header files, for example BLAS.h,
> are often not compatible with gfortran function declarations.  To take
> one small example, in src/main/print.c, we have
>
> void NORET F77_NAME(xerbla)(const char *srname, int *info)
>
> so xerbla_ is defined with two arguments.
>
> However, gfortran passes string lengths as hidden arguments.
> You can see this by compiling the small example
>
> $ cat xer.f
>       SUBROUTINE FOO
>       INTEGER INFO
>       CALL XERBLA ('FOO', INFO)
>       END
> $ gfortran -c -fdump-tree-original xer.f
> $ cat xer.f.004t.original
> foo ()
> {
>   integer(kind=4) info;
>
>   xerbla (&"FOO"[1]{lb: 1 sz: 1}, &info, 3);
> }
>
> so here we have three arguments. This mismatch is flagged
> by -Wlto-type-mismatch, which, for example, yields
>
> print.c:1120:12: note: type 'void' should match type 'long int'
> ../../src/extra/blas/blas.f:357:20: warning: type of 'xerbla' does not
> match original declaration [-Wlto-type-mismatch]
>   357 |          CALL XERBLA( 'DGBMV ', INFO )
>
>
> So, why can gcc's r268992 / r269349 matter? Before these patches,
> gfortran used the variadic calling convention for calling procedures
> outside the current file, and the non-variadic calling convention for
> calling procedures found in the current file.
>
> Because the procedures were all compiled as non-variadic, the caller and
> the calle's signature did not match if they were not in the same
> source file, which is an ABI violation.
>
> This violation manifested itself in https://gcc.gnu.org/PR87689 ,
> where the the problem resulted in crashes on a primary gcc platform,
> POWER.
>
> How can this potentially affect R?  After the fix for PR87689,
> gfortran's calls to external procedures are no longer variadic. It is
> quite possible that, while this "works" most of the time, there
> is a problem with a particular LAPACK routine, the call sequence
> leading up to it or the procedures it calls.
>
> How to fix this problem?  The only clear way I see is to fix this
> on the R side, by adding the string lengths to the prototypes.
> These are size_t (64 bit on 64-bit systems, 32 bit on 32-bit
> systems).  You should then try to make --enable-lto pass
> without any warnings.
>
> Regarding LAPACK itself, the default build system for R builds
> it as a shared library.  Offhand, I did not see any way to
> build a *.a file instead, so I could not use LTO to check
> for mismatched prototypes between R and LAPACK.
>
> Of course, I cannot be sure that this is really the root cause
> of the problem you are seeing,but it does seem to fit quite well.
> I hope this analysis helps in resolving this.
>
> Regards
>
>     Thomas
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

[Thomas König]

Hi Tomas,

thanks a lot for your analysis.  I have created
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=90329
for this, and put you in CC (if your e-mail address
for GCC bugzilla is still current).

Regards

        Thomas

[Berend Hasselman]

Hi, Thomas, Tomas,

Doesn't this issue demonstrate the warning and advice given in the  last paragraph of section 6.6
of the "Writing R Extensions" manual:

<ref>
Passing character strings from C to Fortran or vice versa is not portable (and to Fortran 90 or later is even less so). We have found that it helps to ensure that a C string to be passed is followed by several nuls (and not just the one needed as a C terminator). But for maximal portability character strings in Fortran should be avoided.
</ref>

Avoid.

Berend

> On 3 May 2019, at 19:25, Thomas König <[hidden email]> wrote:
>
> Hi Tomas,
>
> thanks a lot for your analysis.  I have created
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=90329
> for this, and put you in CC (if your e-mail address
> for GCC bugzilla is still current).
>
> Regards
>
> Thomas
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

[peter dalgaard]

The point is that LAPACK uses characters as control arguments in multiple places and we don't write the LAPACK Fortran routines. It has long been known that general character strings was a portability issue but many (not just R people) have thought that length-one character were safe to pass as char* pointers. So "avoid" is not really an option if we want to use LAPACK functionality at all.

Workarounds/solutions include:

- disable certain optimizations -- works for now, but doesn't remove the root cause so seems generally fragile

- "onion-skin" all LAPACK routines to call via a Fortran routine that converts integer arguments to the required character -- possible, but it adds overhead and there are hundreds of routines (and it would be kind of ugly!).

- modify LAPACK itself similarly -- requires naming change of routines as per the license, and there are still hundreds of routines; avoids overhead, but creates maintenance nightmare to keep up with changes in LAPACK

- change all prototypes and calls to follow gfortran calling conventions -- still a lot of work since each char* arguments need to be supplemented by a length going at the end of the arglist. If gfortran was the only compiler around, I'd say this would be the least painful route, but still no fun since it requires changes to a lot of user code (in packages too). It is not clear if this approach works with other Fortrans.

- figure out Fortran2003 specification for C/Fortran interoperability -- this _sounds_ like the right solution, but I don't think many understand how to use it and what is implied (in particular, will it require making changes to LAPACK itself?)

- move towards the LAPACKE C interface -- but that also adds onionskin overhead and ultimately calls Fortran in essentially the same way as R does, so doesn't really solve anything at all (unless we can shift responsibility for sorting things out onto the LAPACK team, but I kind of expect that they do not want it.)

- twist the arms of the gfortran team to do something that keeps old code working. Compiler engineers understandably hate that sort of thing, but I seem to recall some precedent (pointer alignment, back in the dark ages?).

-pd

> On 4 May 2019, at 16:49 , Berend Hasselman <[hidden email]> wrote:
>
> Hi, Thomas, Tomas,
>
> Doesn't this issue demonstrate the warning and advice given in the  last paragraph of section 6.6
> of the "Writing R Extensions" manual:
>
> <ref>
> Passing character strings from C to Fortran or vice versa is not portable (and to Fortran 90 or later is even less so). We have found that it helps to ensure that a C string to be passed is followed by several nuls (and not just the one needed as a C terminator). But for maximal portability character strings in Fortran should be avoided.
> </ref>
>
> Avoid.
>
> Berend
>
>> On 3 May 2019, at 19:25, Thomas König <[hidden email]> wrote:
>>
>> Hi Tomas,
>>
>> thanks a lot for your analysis.  I have created
>> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=90329
>> for this, and put you in CC (if your e-mail address
>> for GCC bugzilla is still current).
>>
>> Regards
>>
>> Thomas
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: [phone]
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

[Thomas König]

Hi Peter,

we (the gfortran team) are currently discussing this at
https://gcc.gnu.org/bugzilla/show_bug.cgi?id=90329 . I
invite everybody who has an interest in this topic to
take part in the discussion there.

> Workarounds/solutions include:
>
> - disable certain optimizations -- works for now, but doesn't remove the root cause so seems generally fragile

That looks like a short-term solutuion that could work (at least for
x86_64 using the standard Unix ABI). And yes, it is fragile.

And whatever other solution people come up with, it will still be
fragile unless the caller and the callee agree.

The root cause is that the Fortran LAPACK routines are called from C
via an incompatible call signature.

> - "onion-skin" all LAPACK routines to call via a Fortran routine that converts integer arguments to the required character -- possible, but it adds overhead and there are hundreds of routines (and it would be kind of ugly!).

I agree.

> - modify LAPACK itself similarly -- requires naming change of routines as per the license, and there are still hundreds of routines; avoids overhead, but creates maintenance nightmare to keep up with changes in LAPACK

I agree that this is not a preferred option.

> - change all prototypes and calls to follow gfortran calling conventions -- still a lot of work since each char* arguments need to be supplemented by a length going at the end of the arglist. If gfortran was the only compiler around, I'd say this would be the least painful route, but still no fun since it requires changes to a lot of user code (in packages too). It is not clear if this approach works with other Fortrans.

The interesting thing is that this convention goes back to at least f2c,
which was modeled on the very first Unix compiler.

> - figure out Fortran2003 specification for C/Fortran interoperability -- this _sounds_ like the right solution, but I don't think many understand how to use it and what is implied (in particular, will it require making changes to LAPACK itself?)

That would actually be fairly easy.  If you declare the subroutines
BIND(C), as in

       subroutine foo(a,b) BIND(C,name="foo_")
       real a
       character*1 b
       end

you will get the calling signature that you already have in your C
sources.

This also has the advantage of being standards compliant, and would be
probably be the preferred method.

> - move towards the LAPACKE C interface -- but that also adds onionskin overhead and ultimately calls Fortran in essentially the same way as R does, so doesn't really solve anything at all (unless we can shift responsibility for sorting things out onto the LAPACK team, but I kind of expect that they do not want it.)

I suspect that they will hit the issue, too.

> - twist the arms of the gfortran team to do something that keeps old code working. Compiler engineers understandably hate that sort of thing, but I seem to recall some precedent (pointer alignment, back in the dark ages?).

We're willing to do reasonable things :-) but so far all of the options
we have come up with have very serious drawbacks (see the link to the
PR at the top). If you come up with a suggestion, we'd be more than
happy to look at it.

I think the best option would really be to use BIND(C).

Regards

        Thomas

[Steve Kargl]

On Sat, May 04, 2019 at 06:05:48PM +0200, peter dalgaard wrote:
>
> - figure out Fortran2003 specification for C/Fortran interoperability
> -- this _sounds_ like the right solution, but I don't think many
> understand how to use it and what is implied (in particular, will it
> require making changes to LAPACK itself?)

This is probably the best solution as it should allow portability
to any Fortran processor that supports F2003 or newer standard.
See the gtk-fortran project.  It has a python program that was
used to generate the needed ISO C interfaces.

--
Steve

[Steve Kargl]

On Sat, May 04, 2019 at 06:42:47PM +0200, Thomas König wrote:

>
> > - figure out Fortran2003 specification for C/Fortran interoperability
> > -- this _sounds_ like the right solution, but I don't think many
> > understand how to use it and what is implied (in particular, will
> > it require making changes to LAPACK itself?)
>
> That would actually be fairly easy.  If you declare the subroutines
> BIND(C), as in
>
>        subroutine foo(a,b) BIND(C,name="foo_")
>        real a
>        character*1 b
>        end
>
> you will get the calling signature that you already have in your C
> sources.
>
> This also has the advantage of being standards compliant, and would be
> probably be the preferred method.
>

With the caveat that one may need to use the VALUE attribute to
account for pass-by-value vs pass-by-reference.

--
Steve

[Thomas König]

Hi Steve,

> With the caveat that one may need to use the VALUE attribute to
> account for pass-by-value vs pass-by-reference.

LAPACK should be all pass by reference, it is old F77-style
code (except that the odd ALLOCATABLE array has snuck in
in the testing routines).

[Tomas Kalibera-2]

On 5/4/19 6:49 PM, Steve Kargl wrote:

> On Sat, May 04, 2019 at 06:42:47PM +0200, Thomas König wrote:
>>> - figure out Fortran2003 specification for C/Fortran interoperability
>>> -- this _sounds_ like the right solution, but I don't think many
>>> understand how to use it and what is implied (in particular, will
>>> it require making changes to LAPACK itself?)
>> That would actually be fairly easy.  If you declare the subroutines
>> BIND(C), as in
>>
>>         subroutine foo(a,b) BIND(C,name="foo_")
>>         real a
>>         character*1 b
>>         end
>>
>> you will get the calling signature that you already have in your C
>> sources.
>>
>> This also has the advantage of being standards compliant, and would be
>> probably be the preferred method.
>>
> With the caveat that one may need to use the VALUE attribute to
> account for pass-by-value vs pass-by-reference.

This seems clean solution, but as I said before not easy, because
currently the tradition is to call the Fortran interface directly from C
(not via any C wrappers). This means one could not substitute
LAPACK/BLAS at dynamic linking time, unless all LAPACK/BLAS
implementations agreed on such a C interface. Now the substitution is
based on the original Fortran interface.

In case of R, if we only used the included reference BLAS/LAPACK, we
could do this, define our wrappers, say "c_dgemm" for "dgemm", change R
to call via that interface, ask maintainers of all packages to change
their code to call via their interface, and this should work with all
Fortran 2003 compilers.

But, R is often used also with optimized BLAS/LAPACK implementations
that can be substituted at dynamic linking time. And there we could do
nothing at R level to help: we cannot generate such wrappers for an
existing LAPACK/BLAS implementation (we don't have the source code, the
compiler, etc).

It would be certainly a good thing if BLAS/LAPACK, with all
implementation and uses, switched to a way that is compliant with
current Fortran standard. But this should best start with the reference
BLAS/LAPACK, continue with other BLAS/LAPACK implementations, and then
with systems using those libraries, including R and its packages.
Unless/before this happens, it would really be great if we could still
use gfortran to build and use this fundamental software library.

Best
Tomas

[Janne Blomqvist-3]

On Mon, May 6, 2019 at 11:55 AM Tomas Kalibera <[hidden email]> wrote:

>
> On 5/4/19 6:49 PM, Steve Kargl wrote:
> > On Sat, May 04, 2019 at 06:42:47PM +0200, Thomas König wrote:
> >>> - figure out Fortran2003 specification for C/Fortran interoperability
> >>> -- this _sounds_ like the right solution, but I don't think many
> >>> understand how to use it and what is implied (in particular, will
> >>> it require making changes to LAPACK itself?)
> >> That would actually be fairly easy.  If you declare the subroutines
> >> BIND(C), as in
> >>
> >>         subroutine foo(a,b) BIND(C,name="foo_")
> >>         real a
> >>         character*1 b
> >>         end
> >>
> >> you will get the calling signature that you already have in your C
> >> sources.
> >>
> >> This also has the advantage of being standards compliant, and would be
> >> probably be the preferred method.
> >>
> > With the caveat that one may need to use the VALUE attribute to
> > account for pass-by-value vs pass-by-reference.
>
> This seems clean solution, but as I said before not easy, because
> currently the tradition is to call the Fortran interface directly from C
> (not via any C wrappers). This means one could not substitute
> LAPACK/BLAS at dynamic linking time, unless all LAPACK/BLAS
> implementations agreed on such a C interface. Now the substitution is
> based on the original Fortran interface.
>
> In case of R, if we only used the included reference BLAS/LAPACK, we
> could do this, define our wrappers, say "c_dgemm" for "dgemm", change R
> to call via that interface, ask maintainers of all packages to change
> their code to call via their interface, and this should work with all
> Fortran 2003 compilers.
>
> But, R is often used also with optimized BLAS/LAPACK implementations
> that can be substituted at dynamic linking time. And there we could do
> nothing at R level to help: we cannot generate such wrappers for an
> existing LAPACK/BLAS implementation (we don't have the source code, the
> compiler, etc).
>
> It would be certainly a good thing if BLAS/LAPACK, with all
> implementation and uses, switched to a way that is compliant with
> current Fortran standard. But this should best start with the reference
> BLAS/LAPACK, continue with other BLAS/LAPACK implementations, and then
> with systems using those libraries, including R and its packages.
> Unless/before this happens, it would really be great if we could still
> use gfortran to build and use this fundamental software library.

Hi,

I don't think modifying your own builtin LAPACK makes sense, as you
mention that would make R incompatible with another LAPACK provided by
the system. And modifying LAPACK upstream by adding BIND(C) wouldn't
work either, as that would break all the existing Fortran code calling
LAPACK (as LAPACK is F77-style implicit interfaces, the Fortran caller
has no knowledge of the interface and thus it must match the compiler
default Fortran ABI).

So the remaining place where this could be fixed would be in your C
prototypes for LAPACK functions, so that they match what LAPACK
expects. Arguably that's the correct approach anyway. I suppose for
this task extending the GFortran -fc-prototypes option to generate C
prototypes for external functions as well would help?

AFAICS, this interface mismatch problem affects other Fortran
compilers as well, just that by sheer luck this has worked mostly so
far (that is, other Fortran compilers also expect a hidden string
length argument, with no exception for length==1 strings). But with
increasingly sophisticated interprocedural optimizations such sins can
no longer be forgiven.

--
Janne Blomqvist

[Tomas Kalibera-2]

On 5/6/19 12:57 PM, Janne Blomqvist wrote:

> On Mon, May 6, 2019 at 11:55 AM Tomas Kalibera <[hidden email]> wrote:
>> On 5/4/19 6:49 PM, Steve Kargl wrote:
>>> On Sat, May 04, 2019 at 06:42:47PM +0200, Thomas König wrote:
>>>>> - figure out Fortran2003 specification for C/Fortran interoperability
>>>>> -- this _sounds_ like the right solution, but I don't think many
>>>>> understand how to use it and what is implied (in particular, will
>>>>> it require making changes to LAPACK itself?)
>>>> That would actually be fairly easy.  If you declare the subroutines
>>>> BIND(C), as in
>>>>
>>>>          subroutine foo(a,b) BIND(C,name="foo_")
>>>>          real a
>>>>          character*1 b
>>>>          end
>>>>
>>>> you will get the calling signature that you already have in your C
>>>> sources.
>>>>
>>>> This also has the advantage of being standards compliant, and would be
>>>> probably be the preferred method.
>>>>
>>> With the caveat that one may need to use the VALUE attribute to
>>> account for pass-by-value vs pass-by-reference.
>> This seems clean solution, but as I said before not easy, because
>> currently the tradition is to call the Fortran interface directly from C
>> (not via any C wrappers). This means one could not substitute
>> LAPACK/BLAS at dynamic linking time, unless all LAPACK/BLAS
>> implementations agreed on such a C interface. Now the substitution is
>> based on the original Fortran interface.
>>
>> In case of R, if we only used the included reference BLAS/LAPACK, we
>> could do this, define our wrappers, say "c_dgemm" for "dgemm", change R
>> to call via that interface, ask maintainers of all packages to change
>> their code to call via their interface, and this should work with all
>> Fortran 2003 compilers.
>>
>> But, R is often used also with optimized BLAS/LAPACK implementations
>> that can be substituted at dynamic linking time. And there we could do
>> nothing at R level to help: we cannot generate such wrappers for an
>> existing LAPACK/BLAS implementation (we don't have the source code, the
>> compiler, etc).
>>
>> It would be certainly a good thing if BLAS/LAPACK, with all
>> implementation and uses, switched to a way that is compliant with
>> current Fortran standard. But this should best start with the reference
>> BLAS/LAPACK, continue with other BLAS/LAPACK implementations, and then
>> with systems using those libraries, including R and its packages.
>> Unless/before this happens, it would really be great if we could still
>> use gfortran to build and use this fundamental software library.
> Hi,
>
> I don't think modifying your own builtin LAPACK makes sense, as you
> mention that would make R incompatible with another LAPACK provided by
> the system. And modifying LAPACK upstream by adding BIND(C) wouldn't
> work either, as that would break all the existing Fortran code calling
> LAPACK (as LAPACK is F77-style implicit interfaces, the Fortran caller
> has no knowledge of the interface and thus it must match the compiler
> default Fortran ABI).

I meant creating another layer of functions using BIND(C), which will
have different names (e.g. c_dgemm), and call into the original Fortran
functions (e.g. dgemm). Fortran programs could still call into the
original Fortran functions. Something like c_dgemm in Chapter 10.2 of
"Numerical Computing with Modern Fortran" by Richard J. Hanson and Tim
Hopkins, the code is available from
https://archive.siam.org/books/ot134/Chapter10/index.php. So
LAPACK/BLACK implementations won't break their applications by adding
these new functions. Also indeed it would be useful to have C prototypes
matching those interfaces distributed with BLAS/LAPACK.
> So the remaining place where this could be fixed would be in your C
> prototypes for LAPACK functions, so that they match what LAPACK
> expects. Arguably that's the correct approach anyway. I suppose for
> this task extending the GFortran -fc-prototypes option to generate C
> prototypes for external functions as well would help?

Well but the non-BIND(C) interface is different for different Fortran
compilers, and it changed even between gfortran versions (the type of
the lengths). So we would not be able to switch BLAS/LAPACK
implementations anymore at dynamic linking time. We would still have to
talk to all R package maintainers that manage packages that call
directly to BLAS/LAPACK.

I still think -fc-prototypes would be useful, certainly it should work
for all BIND(C) procedures (the 20190426 does not for the dgemm example
from the book), and it would be useful even for non-BIND(C) procedures.

> AFAICS, this interface mismatch problem affects other Fortran
> compilers as well, just that by sheer luck this has worked mostly so
> far (that is, other Fortran compilers also expect a hidden string
> length argument, with no exception for length==1 strings). But with
> increasingly sophisticated interprocedural optimizations such sins can
> no longer be forgiven.

Best
Tomas

[Thomas Koenig-6]

Hi,

gfortran trunk and 9-branch now have an option to automatically
generate C prototypes for old-style F77 procedures.  I just did

for a in *.f; do gfortran -fsyntax-only -fc-prototypes-external $a >
${a%.f}.h; done

in the src/modules/lapack directory.  This generates header
files which contain prototypes like

int ilaenv_ (int *ispec, char *name, char *opts, int *n1, int *n2, int
*n3, int *n4, size_t name_len, size_t opts_len);
void dlacn2_ (int *n, double *v, double *x, int *isgn, double *est, int
*kase, int *isave);
void dlaln2_ (int_least32_t *ltrans, int *na, int *nw, double *smin,
double *ca, double *a, int *lda, double *d1, double *d2, double *b, int
*ldb, double *wr, double *wi, double *x, int *ldx, double *scale, double
*xnorm, int *info);
void dlabad_ (double *small, double *large);
void drscl_ (int *n, double *sa, double *sx, int *incx);
void dlatrs_ (char *uplo, char *trans, char *diag, char *normin, int *n,
double *a, int *lda, double *x, double *scale, double *cnorm, int *info,
size_t uplo_len, size_t trans_len, size_t diag_len, size_t normin_len);

which could serve as the basis for adjusting the calling sequence
for the C bindings to what the compiler expects.

I checked, and it appears that at least ifort uses the same convention
as gfortran 8/9 regarding character argument passing.

Regards

        Thomas

[Steve Kargl]

On Sat, May 11, 2019 at 05:04:06PM +0200, Thomas Koenig wrote:

>
> gfortran trunk and 9-branch now have an option to automatically
> generate C prototypes for old-style F77 procedures.  I just did
>
> for a in *.f; do gfortran -fsyntax-only -fc-prototypes-external $a >
> ${a%.f}.h; done
>
> in the src/modules/lapack directory.  This generates header
> files which contain prototypes like
>
> int ilaenv_ (int *ispec, char *name, char *opts, int *n1, int *n2, int
> *n3, int *n4, size_t name_len, size_t opts_len);

Any chance that this can be created without the parameter names?

For example,

int ilaenv_ (int *, char *, char *, int *, int *, int *, int *,
        size_t , size_t );

For bonus points, a tab on continuation lines would be nice.

--
Steve

[Thomas Koenig-6]

Hi Steve,

> On Sat, May 11, 2019 at 05:04:06PM +0200, Thomas Koenig wrote:
>>
>> gfortran trunk and 9-branch now have an option to automatically
>> generate C prototypes for old-style F77 procedures.  I just did
>>
>> for a in *.f; do gfortran -fsyntax-only -fc-prototypes-external $a >
>> ${a%.f}.h; done
>>
>> in the src/modules/lapack directory.  This generates header
>> files which contain prototypes like
>>
>> int ilaenv_ (int *ispec, char *name, char *opts, int *n1, int *n2, int
>> *n3, int *n4, size_t name_len, size_t opts_len);
>
> Any chance that this can be created without the parameter names?
>
> For example,
>
> int ilaenv_ (int *, char *, char *, int *, int *, int *, int *,
> size_t , size_t );

I never understood why people leave out the variable names in
prototypes.

Having said that, I think we should leave them in for the BIND(C)
routines, because people (like me :-) will use them for writing
their own C functions, then it is good to have the argument
names there from the start.

For the external procedures, sure, we can make them optional.

> For bonus points, a tab on continuation lines would be nice.

Currently, there is no line break in the generated prototypes,
so also no tab.  So yes, this probably be a good idea.

Regards

        Thomas