I attempted to reproduce the BLAS+F2C bug on Windows 11-64, with some success. I could not use RonShepard’s F90 source as input to F2C, and I did not want to cloud the question by using F77 code with incorrect type declarations for BLAS routine arguments or BLAS function return values. Here is my F77 test code:
program test_blas
implicit none
real dot, snorm, sdot, snrm2, x(3), y(3)
character*20 cfmt
data x/1.0, 2.0, 3.0/, y /3.0, 0.0, 4.0/
data cfmt/'(2(a6,f10.3))'/
dot = sdot(3,x,1,y,1)
snorm = snrm2(3,y,1)
write(*,cfmt) ' sdot=', dot, ' snrm2=', snorm
end
Just now, I downloaded the current F2C sources from Netlib and built f2c.exe and vcf2c.lib using “Microsoft (R) C/C++ Optimizing Compiler Version 19.33.31629 for x64”. The above F77 source, converted using the freshly built f2c.exe, gave this C source:
/* rssp.f -- translated by f2c (version 20200916).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Table of constant values */
static integer c__3 = 3;
static integer c__1 = 1;
/* Main program */ int MAIN__(void)
{
/* Initialized data */
static real x[3] = { 1.f,2.f,3.f };
static real y[3] = { 3.f,0.f,4.f };
static char cfmt[20] = "(2(a6,f10.3)) ";
/* Builtin functions */
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
/* Local variables */
static real dot;
extern doublereal sdot_(integer *, real *, integer *, real *, integer *),
snrm2_(integer *, real *, integer *);
static real snorm;
/* Fortran I/O blocks */
static cilist io___6 = { 0, 6, 0, cfmt, 0 };
dot = sdot_(&c__3, x, &c__1, y, &c__1);
snorm = snrm2_(&c__3, y, &c__1);
s_wsfe(&io___6);
do_fio(&c__1, " sdot=", (ftnlen)6);
do_fio(&c__1, (char *)&dot, (ftnlen)sizeof(real));
do_fio(&c__1, " snrm2=", (ftnlen)7);
do_fio(&c__1, (char *)&snorm, (ftnlen)sizeof(real));
e_wsfe();
return 0;
} /* MAIN__ */
/* Main program alias */ int test_blas__ () { MAIN__ (); return 0; }
The bug is brazen – the line with doublereal. With this bug left unfixed, the output from the C code is wrong:
sdot= .000 snrm2 25.000
Replacing “doublereal” by “real” in the F2C-generated C file fixes the bug, and the output becomes
sdot= 15.000 snrm2 5.000
In this reconstructed scenario, then, the culprit is F2C. Given that F2C is obsolete, we should not think of fixing the F2C sources to address this problem. Whoever (Apple, in 2022?) uses/used F2C and a C compiler to build Accelerate needs to do the right thing – recompile BLAS, Lapack, etc., from Fortran sources or fix their F2C-converted C sources. Furthermore, at least in this reconstructed scenario, the conclusion “This is a C language error, one that has propagated now for several decades.” does not apply.
P.S. Line 81 of file main.c of the F2C converter sources contains:
int forcedouble = YES; /* force real functions to double */
Changing YES to NO on this line and rebuilding f2c.exe, I found that this change removes the problem with the test code, produces the correct type declarations, and the resulting C code when built produces the correct result. Of course, there is no telling what things break with such a change to a very complicated body of code such as that of F2C.
P.P.S. Many of us do not care to read old manuals carefully. I’m guilty as well. Here are a few lines from the the third page of the F2C document:
A function of type integer, logical, or double precision must be declared as a C function that returns the corresponding type. If the -R option is in effect (see Appendix B), the same is true of a function of type real; otherwise, a real function must be declared as a C function that returns doublereal; this hack facilitates our VAX regression testing, as it duplicates the behavior of our local Fortran compiler (f77).
Thus, to cause F2C to produce a C function code whose return variable is float from Fortran 77 source of a function of type real, one should use the -R option of F2C. Not being aware of that flag is easier to excuse than AT&T’s making the default behavior of F2C to be to deliver an unwanted change of the return type of a single precision real function in Fortran to a double precision function in C.