On Tuesday, September 12, 2023 at 7:36:32 PM UTC+3, db wrote:
> On 12.09.2023 16.32, Michael S wrote:
> > On Tuesday, September 12, 2023 at 4:23:09 PM UTC+3, Ben Bacarisse wrote:
> >> Michael S <already...@yahoo.com> writes:
> >>
> >>> On Tuesday, September 12, 2023 at 2:26:49 AM UTC+3, Ben Bacarisse wrote:
> >>>> db <dieterh...@gmail.com> writes:
> >>>>
> >>>>> On 10.09.2023 15.22, Ben Bacarisse wrote:
> >>>>>> db <dieterh...@gmail.com> writes:
> >>>>>>
> >>>>>>> On 09.09.2023 11.02, Mut...@dastardlyhq.com wrote:
> >>>>>>>> On Sat, 9 Sep 2023 10:54:33 +0200
> >>>>>>>> db <dieterh...@gmail.com> wrote:
> >>>>>>>>> but now it says that C-double doesn't exist.
> >>>>>>>>>
> >>>>>>>>> Any ideas?
> >>>>>>>> C++ mangles function names so that the name you give the function isn't
> >>>>>>>> the
> >>>>>>>> name the dynamic loader sees in the binary. You need to put this around whatever
> >>>>>>>> function you want to call from Fortran so it has C style linkage which means
> >>>>>>>> the name remains the same:
> >>>>>>>> extern "C" {
> >>>>>>>> void calledFromFortran()
> >>>>>>>> {
> >>>>>>>> :
> >>>>>>>> }
> >>>>>>>> }
> >>>>>>>>
> >>>>>>> Where should this go in the code? I added it at the end as is,
> >>>>>>> and now I get this:
> >>>>>>>
> >>>>>>> ~/ownlib90/tests> ./jbandtest
> >>>>>>> jband.cpp:88:1: error: expected initializer before 'extern'
> >>>>>>> 88 | extern "C" {
> >>>>>>> | ^~~~~~
> >>>>>>> jband.cpp:96:1: error: expected unqualified-id before '{' token
> >>>>>>> 96 | {
> >>>>>>> | ^
> >>>>>>> jband.cpp: In function 'void calledFromFortran()':
> >>>>>>> jband.cpp:350:9: error: expected primary-expression before ':' token
> >>>>>>> 350 | :
> >>>>>>> ...
> >>>>>> You are giving people very little to go on. There's an unknown command
> >>>>>> (./jbandtest) the gives some errors where only one line of code can be
> >>>>>> seen.
> >>>>>> Can you produce a cut-down version to the program that shows the same
> >>>>>> error -- mimimal example of the C++ and of the Fortran that should
> >>>>>> calls it but which fails? That would mean you could post the code.
> >>>>>> If that is not possible, can you link to the source so those who know
> >>>>>> both C++ and Fortran can take a proper look?
> >>>>>
> >>>>> I put the C++ function here: www.dieterbritz.dk/jband.cpp
> >>>>> After getting errors about the extern lines, whch I had
> >>>>> added at the end, I shifted it to near the start:
> >>>>>
> >>>>> #include <iostream>
> >>>>> #include <fstream>
> >>>>> #include <iomanip>
> >>>>> #include "math.h"
> >>>>>
> >>>>> using namespace std;
> >>>>>
> >>>>> extern "C" {
> >>>>>
> >>>>> void calledFromFortran()
> >>>>> {
> >>>>> :
> >>>>> }
> >>>>>
> >>>>> }
> >>>> Remove them all. You have two options. One: correct the C++. To do
> >>>> this put this, and only this,
> >>>>
> >>>> extern "C" {
> >>>> const long double iband(const long double tbar);
> >>>> }
> >>>>
> >>>> just after the #include lines. This is perfect (the first const is a
> >>>> meaningless word here and the compiler may warn you about that) but it
> >>>> will do.
> >>>>
> >>>> Next replace
> >>>>
> >>>> #include "math.h"
> >>>>
> >>>> with
> >>>>
> >>>> #include <math.h>
> >>>>
> >>>> I'd also remove the two unnecessary #include lines (for iomanip and
> >>>> fstream). All of these are "tidying up". Only the extern "C" part I
> >>>> gave corrects an actual mistake.
> >>>>
> >>>> The second option is the turn the code in C which links directly with
> >>>> Fortran. The only C++ used by the code is one "cout << " line (not a
> >>>> good way to report an error, but that's not important right now) and a
> >>>> setting for pi that is not a compile-time constant (4*atanl(1)) written
> >>>> in a more fussy way.
> >>>>
> >>>> To convert to C, copy the file to jband.c (rather than .cpp) and replce
> >>>> the declaration for pi. Given the code is littered with explicit 21
> >>>> digits decimals (a little more than the precision of C's usualy long
> >>>> double type) I would just write
> >>>>
> >>>> static const long double pi = 3.14159265358979323846L;
> >>>>
> >>>> Now replace the cout line with
> >>>>
> >>>> fputs("\nError in stredges::iband(): tbar must be > 0", stderr);
> >>>>
> >>>> This writes the message to the error stream and put only these two
> >>>> include lines at the top:
> >>>>
> >>>> #include <math.h>
> >>>> #include <stdio.h>
> >>>>
> >>>> Now the file is C code.
> >>>>
> >>>> Frankly, this is the way I'd go, since linking C is just that bit
> >>>> simpler than linking C++. But if the code is going to become fully
> >>>> fledged (and complex) C++, obviously use my first option.
> >>>>> I still get errors pertaining to these lines. But after adding
> >>>>> the -c to the two compilation lines in the script, these are
> >>>>> now the only errors. I think I am getting closer. The errors
> >>>>> are now
> >>>> Yes, you must compile it like this:
> >>>>
> >>>> g++ -c jband.cpp (for C++ code)
> >>>> gcc -c jband.c (for C code)
> >>>>
> >>>> You don't need any -o but you could add in lots of warnings like
> >>>>
> >>>> gcc -Wall -Wextra -c jband.c
> >>>>> jband.cpp: In function 'void calledFromFortran()':
> >>>>> jband.cpp:12:9: error: expected primary-expression before ':' token
> >>>>> 12 | :
> >>>>> | ^
> >>>>> jband.cpp: At global scope:
> >>>>> jband.cpp:96:1: error: expected initializer before 'extern'
> >>>>> 96 | extern "C" {
> >>>>> | ^~~~~~
> >>>>> jband.cpp:104:1: error: expected unqualified-id before '{' token
> >>>>> 104 | {
> >>>>> | ^
> >>>>> jbandtest.f90:28:14:
> >>>>>
> >>>>> 28 | FORTJBAND = JBAND (tbar)
> >>>>> | 1
> >>>>> Error: Type mismatch in argument 'tbar' at (1); passed REAL(16) to REAL(10)
> >>>>>
> >>>>> The last error will go away together with the others.
> >>>> Michael S has given you one way to go for that. My comments were going
> >>>> to be more general. You have several ways of specifying the types, and
> >>>> I could not tell if you really wanted all of them. In general, I would
> >>>> define the Fortran types from the C ones (since those are the types you
> >>>> need for the call and return) and leave it at that.
> >>>>
> >>>> But the big one is that the Fortran code seems to think the function is
> >>>> called jband, when in fact it's called iband.
> >>>>
> >>>> I went for this:
> >>>>
> >>>> module STUFF
> >>>> use iso_c_binding
> >>>> integer, parameter :: qud=real(c_long_double)
> >>>> end module STUFF
> >>>> program JBAND_TEST
> >>>> use STUFF; implicit none
> >>>> real(qud) :: tbar, curr, FORTJBAND
> >>>> do
> >>>> read *, tbar
> >>>> if (tbar < 0) exit
> >>>> curr = FORTJBAND (tbar)
> >>>> print '(" curr =", f10.4)', curr
> >>>> enddo
> >>>> end program JBAND_TEST
> >>>>
> >>>> function FORTJBAND (tbar)
> >>>> use STUFF; implicit none
> >>>> interface
> >>>> function IBAND (tbar) bind(c)
> >>>> use STUFF
> >>>> real(qud), VALUE :: tbar
> >>>> real(qud) :: IBAND
> >>>> end function iband
> >>>> end interface
> >>>
> >>> I think, you are taking chances here. Even if it happens to work
> >>> on x86-64 Linux, it can fail on some other target where calling
> >>> conventions for long double return value in 'C' differ from those
> >>> of Fortran. The method shown in one of my previous posts
> >>> appear more robust:
> >>>
> >>> interface
> >>> real(c_long_double) function IBAND (tbar) bind(c)
> >>> import c_long_double
> >>> real(c_long_double), value :: tbar
> >>> end function IBAND
> >> I think I see -- keep the C-related types in the smallest interface? I
> >> didn't study your example, but the OP renamed the C function and gave
> >> the resulting Fortran function a native Fortran type. Is that what you
> >> advise?
> >>
> >> The OP should have taken your Fortran and applied my C++ edits. Maybe
> >> they took my Fortran and it did not work on Windows. That would be a
> >> shame.
> >>>> BTW, you said the C++ was written by a C++ expert. I don't think it
> >>>> was!
> >>>
> >>> With no relationships to C, C++ or Fortran, I don't like the way they
> >>> evaluate their polynomials from numerical perspective . But I assume
> >>> that algorithm was copied 'as is' from 40 y.o. book so it's not a
> >>> fault of C++ programmer.
> >> I was just remarking on the C++ itself which has some 'non-expert'
> >> attributes.
> >>> Despite my disliking, their method is good enough to achieve their modest
> >>> specified goal for relative accuracy: of 14-15 significant digits.
> >>> With smarter method the same or better accuracy can be achieved using
> >>> standard 'double" arithmetic (IEEE-754 binary64) thus making the routine
> >>> portable to platforms that have no support for extended precision
> >>> math.
> >> Interesting. I don't know the function being evaluated at all.
> >>
> >
> > I don't know it either. But I had drawn the graphs and looked at
> > coefficients in two ranges where precision is worst: [0.0081: 0.0158] and
> > [0.0158:0.0631].
> > A simple affine transform of the polynomial could either improve precision
> > of result by 4 orders of magnitude or reduce needs for precision of arithmetic
> > by 3 orders of magnitude while still somewhat improving the precision of result.
> > I'd guess, more advanced techniques can improve it further in both dimensions,
> > but low hanging fruits first.
> >
> >> --
> >> Ben.
> The person who gave me the C++ code is an older experienced C++
> programmer. Maybe he uses out of data C++ features but the
> function works as is (without the extra lines given to me here,
> which were meant to help with calling it from Fortran). He also
> gave me a little C++ program to call it, and it produced
> output that I was able to verify.
>
> The problem is electrochemical, the current vs time at a flat band
> electrode flush with an insulating plane, for which there is no
> analytical solution, and the function gives highly accurate values
> using piecewise polynomial fits. Actually, he is over the top with
> the accuracy because in practice, comparing with experiments, 3-4
> figures are all that we need but he likes to do the best that can
> be done.
>
> --
> Dieter Britz

If all you need is 3-4 significant figures then you can use much simpler
calculation. Like that (pay attention, it's C rather than C++).

#include <math.h>
#include <stdio.h>

static double polyval(const double coef[], int order, double x);

double iband(double tbar)
{ if (tbar <= 0) {
printf("\nError in stredges::iband(): tbar must be > 0");
return 0;
}

// Coefficients of the polynomial approximation for tbar < 0.76
static const double c1[] = {
-2.0610398523438721606087449e-01,
+3.8436669168692465117185331e-01,
-2.0618169531912929413418346e-01,
-2.6474292774634270832496992e-02,
+3.1131226184858380642960328e-03,
+5.6415965971931841142037414e-01,
};
const int c1_order = sizeof(c1)/sizeof(c1[0])-1;

if (tbar < 0.76)
return polyval(c1, c1_order, tbar)/sqrt(tbar) + 1;

// Coefficients of the polynomial approximation for tbar >= 0.76
static const double c2[] = {
+1.3247948853085559402034700e+01,
-2.8936433794424133211569633e+01,
+2.2857148655763359956491203e+01,
-5.3743813763778698859386414e+00,
-3.5789894047006239277645653e+00,
+6.2827953597648754628587263e+00,
};
const int c2_order = sizeof(c2)/sizeof(c2[0])-1;
static const double c2_b2 = 3.00445175355660613529;
static const double c2_b1 = 1.93;

double log_t = log(tbar);
return polyval(c2, c2_order, 1/(log_t+c2_b1))/(log_t+c2_b2);
}

static double polyval(const double coef[], int order, double x)
{ double sum = coef[0];
for (int i = 0; i < order; ++i)
sum = sum * x + coef[i+1];
return sum;
}