topic I see! Then I fix my code to in Intel® Fortran Compiler

Some Problems about using ODRPACK95

yan__san-chi — Sat, 31 Dec 2016 08:48:33 GMT

Hi!

I have some questions about passing arguments from C to Fortran with ISO_C_BINDING.

1. In ODRPACK, there are arguments need to pass 3D array, for example in wrapper code with ISO_C_BINDING for ODRPACK95 :

!!!!!Optional variable
    !!!!!!!!Unfinished
    
    integer(c_int), optional :: IFIXB(1:NP),IFIXX(:,:),JOB,NDIGIT,MAXIT&
    ,IPRINT,LUNERR,LUNRPT,INFO
    
    real(c_double),  optional :: WE(1:N,1:NQ,1:NQ),WD(1:N,1:M,1:M),&
        TAUFAC,SSTOL,PARTOL,STPB(1:NP),STPD(:,:),&
        SCLB(1:NP),SCLD(:,:),LOWER(1:NP),UPPER(1:NP)
    
     integer(c_int), optional,pointer :: IWORK(:)
    real(c_double), optional :: DELTA(1:N,1:M),WORK(:)

In C, how to pass WD to Fortran? My C code shown below

extern "C" {

    void wrapper_ODR(void(*f)(int *N, int *M, int *NP, int *NQ, int *LDN, int *LDM, int *LDNP, \
        double BETA[2], double XPLUSD[4][1], int IFIXB[2], int IFIXX[1][4], int *LDIFX, \
        int *IDEVAL, double F[1][4], double FJACB[1][2][4], double FJACD[1][1][4], int*ISTOP), \

        int N,int M,int NP,int NQ,double BETA[],double Y[][4],double X[][4],\
        double DELTA[][4],double WE[],double WD[1][1][4],int IFIXB[],int IFIXX[],\
        int *JOB,int *NDIGIT,double *TAUFAC,double *SSTOL, double *PARTOL,\
        int *MAXIT, int *IPRINT, int *LUNERR, int *LUNRPT,double STPB[],\
        double STPD[], double SCLB[], double SCLD[], double WORK[], double IWORK[],\
        int *INFO, double LOWER[2], double UPPER[]);

}


int main(){

	int NP = 2, N = 4, M = 1, NQ = 1;
	double BETA[2] = { 2.0, 0 };
	double X[1][4] = { 0.0, 1.0, 2.0, 3.0 };
	double Y[1][4] = { 2.0, 5.0, 8.0, 11.0 };
	double UPPER[2] = { 2.5, 1 };
	double WD[1][1][4];

	for (int i = 0; i < 4; i++){
		if (i == 0){
			WD[1][1] = 0.05;
		}
		else {

			WD[1][1] = 1;
		}
	}

	wrapper_ODR(&FCN, N, M, NP, NQ, BETA, Y, X, NULL, NULL, WD, NULL, NULL, NULL, NULL, NULL, \
		NULL, NULL,NULL, NULL, NULL, NULL, NULL,NULL, NULL, NULL, NULL, NULL, NULL, \
		NULL, UPPER);

In extern "C" function prototype, WD is declared 3D array with size, and in "int main",

I set value for WD then call wrapper_ODR. But results(as attached file) show it's not correct I expected.

It looks like I have wrong code for 3D array, well actually I don't know how to pass 3D array to Fortran.

2. Data type of argument "DELTA" is array with pointer in ODRPACK, how to use it in C?

Unlike the older ODRPACK, the

mecej4 — Sat, 31 Dec 2016 17:11:00 GMT

Unlike the older ODRPACK, the main ODRPACK-95 subroutine takes quite a few arguments that are arrays of assumed shape. To call ODRPACK-95 from C, you would need to pass not just the base address of the array and leading dimension(s), but descriptors.

This can certainly be done if necessary, but doing so correctly needs that you know the structure of these descriptors, use the header file ISO_Fortran_binding.h containing declarations of the descriptor structures, and assign correct values to all the components of all the descriptors. This is considerably more complicated than to call ODRPACK-95 from Fortran or to call ODRPACK-77 from C.

The function argument, FCN, however, is more straightforward to write in C, since it does not have any assumed-shape arguments. Here, for example, is the C version of FCN for the DRIVE1.F example in the ODRPACK-95 distribution.

#include <math.h>

void FCN(int *pn, int *pm, int *pnp, int *pnq, int *pldn, int *pldm, int *pldnp,
         double *pbeta, double *pxplusd,
         int *pifixb, int *pifixx, int *pldifx,
         int *pideval, double *pf, double *pfjacb,double *pfjacd,
         int *istop){
int i, l, nq=*pnq, n=*pn, ldn=*pldn, ldm=*pldm, ldnp=*pldnp, ideval=*pideval;
double *beta=pbeta-1;
if(beta[1] < 0.0){
        *istop = 1; return;
        }
else istop = 0;
if(ideval % 10)
        for(l=0; l<nq; l++)for(i=0; i<n; i++){
           double xp = exp(beta[3]*pxplusd) - 1.0;
           pf[i+l*ldn] = beta[1] + beta[2]*xp*xp;
           }
if((ideval/10) % 10)
        for(l=0; l<nq; l++)for(i=0; i<n; i++){
           double xp = exp(beta[3]*pxplusd);
           pfjacb[i+l*ldn*ldnp] = 1.0;
           pfjacb[i+ldn+l*ldn*ldnp] = (xp-1.0)*(xp-1.0);
           pfjacb[i+2*ldn+l*ldn*ldnp] = 2*beta[2]*(xp-1.0)*xp*pxplusd;
           }
if((ideval/100) % 10)
        for(l=0; l<nq; l++)for(i=0; i<n; i++){
           double xp = exp(beta[3]*pxplusd);
           pfjacd[i+l*ldn*ldm] = 2*beta[2]*(xp-1.0)*xp*beta[3];
           }
return;
}

Hi!

yan__san-chi — Sun, 01 Jan 2017 10:33:23 GMT

Hi!

2D array like X,Y I can pass C array to Fortran, but in the same way for passing 3D array WD to Fortran above, I got unusual(as attached file) value. I don't know why?

Quote:yan, san-chi wrote: 2D

mecej4 — Sun, 01 Jan 2017 17:47:13 GMT

yan, san-chi wrote:
2D array like X,Y I can pass C array to Fortran, but in the same way for passing 3D array WD to Fortran above, I got unusual(as attached file) value. I don't know why?

I already stated in #2 why this will not work. The main file odr.f of ODRPACK-95 applies the intrinsic function SIZE to several assumed shape array arguments in many places. Since your C program passes only the base address of the array, instead of a descriptor as the Fortran subroutine expects, the latter will pull a wrong value for the size of the array and compute garbage results after that.

I see! Then I fix my code to

yan__san-chi — Tue, 03 Jan 2017 12:36:15 GMT

I see! Then I fix my code to pass 3D array. It works! Thanks for teaching!