CODEWITHSUNDEEP
c++arrayspointersparametersfortran
user1666287
user1666287

Reputation: 11

Passing C++ function pointers to Fortran subroutine

I have been trying to debug this code which calls a Fortran subroutine, DLSODE. Two of the parameters of this subroutine are functions, FEX and JEX, which I am passing as function pointers. From what I understand, any of the functions used by a fortran subroutine must have call by reference parameters, as fortran does not accept call by value.

My problem is the error I get seems to indicate I have a discrepancy between the data types I declared in the function declaration of the pointer and the function definition below, and I have not been successful in debugging it.

#include <iostream>
#include <cmath>
using namespace std;




extern"C" {void DLSODE_(void (*fex)(int, double, double[], double[]), int *NEQ,
                        double *Y[], double *T, double *TOUT, int *ITOL,
                        double *RTOL, double *ATOL[], int *ITASK, int *ISTATE,
                        int *IOPT, double *RWORK[], int *LRW, int *IWORK[],
                        int *LIW, void (*jex)(int, double, double[], int, 
                                              int, double [3][3], int), 
                        int *MF);
}    




void FEX (int &NEQ, double &T, double Y[], int YDOT[]);




void JEX (int &NEQ, double &T, double Y[], int &ML, int &MU, double PD[3][3], 
          int &NRPD);




int main(){

    int IOPT, IOUT, ISTATE, ITASK, ITOL, IWORK[23], LIW, LRW, MF, NEQ, ICOUNT=1;
    double ATOL[3], RTOL, RWORK[58], T, TOUT, Y[3];
    NEQ = 3;
    Y[0] = 1;
    Y[1] = 0;
    Y[2] = 0;
    T = 0;
    TOUT = 0.4;
    ITOL = 2;
    RTOL = pow(10,-4);
    ATOL[0] = pow(10,-6);
    ATOL[1] = pow(10,-10);
    ATOL[2] = pow(10,-6);
    ITASK = 1;
    ISTATE = 1;
    IOPT = 0;
    LRW = 58;
    LIW = 23;
    MF = 21;

    for (ICOUNT; ICOUNT <13; ICOUNT++)
    {
        DLSODE_(&FEX, &NEQ, &Y, &T, &TOUT, &ITOL, &RTOL, &ATOL, &ITASK, &ISTATE, 
                &IOPT, &RWORK, &LRW, &IWORK, &LIW, &JEX, &MF);
        cout<<"At t= "<<T<<" y= "<<Y[1]<<" "<<Y[2]<<" "<<Y[3]<<"\n";
    }
}




void FEX (int &NEQ, double &T, double Y[], double YDOT[]){
    YDOT[0] = -.04*Y[0] + pow(10,4)*Y[1]*Y[2];
    YDOT[2] = 3*pow(10,7)*Y[1]*Y[1];
    YDOT[1] = -YDOT[0] - YDOT[2];
}




void JEX (int &NEQ, double &T, double Y[], int &ML, int &MU, double &PD[3][3], 
          int &NRPD){
    PD[0][0] = -.04;
    PD[0][1] = 1*pow(10,4)*Y[2];
    PD[0][2] = 1*pow(10,4)*Y[1];
    PD[1][0] = .04;
    PD[1][2] = -PD[0][2];
    PD[2][1] = 6*pow(10,7)*Y[1];
    PD[1][1] = -PD[0][1] - PD[2][1];
}

An the error message is

 cversion1.cpp: In function ‘int main()’:
 cversion1.cpp:44:52: error: invalid conversion from ‘void (*)(int&, double&, double*, int*)’ to ‘void (*)(int, double, double*, double*)’
 cversion1.cpp:44:52: error: cannot convert ‘double (*)[3]’ to ‘double**’ for argument ‘3’ to ‘void DLSODE_(void (*)(int, double, double*, double*), int*, double**, double*, double*, int*, double*, double**, int*, int*, int*, double**, int*, int**, int*, void (*)(int, double, double*, int, int, double (*)[3], int), int*)’
 cversion1.cpp: At global scope:
 cversion1.cpp:55:77: error: declaration of ‘PD’ as array of references
 cversion1.cpp:55:78: error: expected ‘)’ before ‘,’ token
 cversion1.cpp:56:4: error: expected unqualified-id before ‘int’

Thank you,

Upvotes: 1

Views: 1538

Answers (2)

M. S. B.
M. S. B.

Reputation: 29401

The modern way to interface C or C++ to Fortran is to use the Fortran ISO C Binding. This provides way to interface Fortran and C (and C++ via 'extern "C"') that is part of the Fortran language standard and is therefore compiler and platform independent. Statements such as "Fortran passes by reference" and that Fortran compilers append underscores (made in in the webpage linked in the other answer) are no longer necessary. The programmer can specify calling conventions compatible with C and the precise names. This is if you are allowed to modify the Fortran code. Or you can write a glue routine in Fortran between the C/C++ routine and the original Fortran routine. There are many previous questions/answers on Stackoverflow about the ISO C Binding. There are examples in the gfortran manual.

Upvotes: 4

NovaDenizen
NovaDenizen

Reputation: 5325

void (*fex)(int, double, double[], double[])

does not match 

void FEX (int &NEQ, double &T, double Y[], double YDOT[])

You've declared argument 3 of DLSODE to be double *Y[], an array of pointers to double, but this should just be double Y[]. You're passing in &Y, and Y is declared as double Y[3], but you should just pass in Y.

In your JEX declaration, you have double &PD[3][3], but you just want double PD[3][3] here.

You should express calculations like 3*pow(10,7) as 30000000 or 3E7, which are compile-time constants.

Also keep in mind that C and C++ use row-major order for 2-d arrays, but FORTRAN uses column-major order. This will have the effect of transposing arrays when you pass them across the C/FORTRAN boundary.

That covers the error messages you've posted and all the errors I've seen.

Also, this seems like a good reference for what you're trying to do.

Upvotes: 3

Related Questions