Wrong values with Fortran subroutine called from C++

Question

I'm trying to call a Fortran subroutine from C++. As a reference I have a test program in Fortran calling the same subroutine, which gives correct values for the procedure. I have extremely limited Fortran knowledge, so I can't figure out why am I obtaining garbage values.

The Fortran code is as follows:

program test_9j

  implicit none

  integer :: ll1, ll2, ll3, ll4, ll5, ll6, ll7, ll8, llmin, llmax, ndim, ier, i

  real(8) :: L1, L2, L3, L4, L5, L6, L7, L8, LMIN, LMAX

  real(8), allocatable, dimension(:) :: ninecof

  ll1 = 2
  ll2 = 2
  ll3 = 2
  ll4 = 2
  ll5 = 2
  ll6 = 2
  ll7 = 2
  ll8 = 2

  llmin = 0
  llmax = 6

  print*, llmin, llmax, llmax-llmin+1
  print*, ""

  L1 = real(ll1,kind=8)
  L2 = real(ll2,kind=8)
  L3 = real(ll3,kind=8)
  L4 = real(ll4,kind=8)
  L5 = real(ll5,kind=8)
  L6 = real(ll6,kind=8)
  L7 = real(ll7,kind=8)
  L8 = real(ll8,kind=8)

  LMIN = real(llmin,kind=8)
  LMAX = real(llmax,kind=8)

  ndim = llmax - llmin + 1

  ALLOCATE(ninecof(ndim),stat=ier)

!  print*, ier

  ninecof = 0.0

  print*, "main Hi1"

!  print*, "LL", LMIN, LMAX
  CALL w9j(ll1,ll2,ll3,ll4,ll5,ll6,ll7,ll8,llmin,llmax,ndim,ninecof)

!  print*, "main Hi2"

  do i=1,ndim
    print*, ninecof(i)
  end do

  print*, "main Hi3a"

  !DEALLOCATE(ninecof)

  print*, "main Hi3b"

end program test_9j

This gives a series of correct outputs for the subroutine. When I try this in C++ however, it doesn't work: (note: w9j subroutine calls drc6j subroutine which is why I included it in the extern definition. drc6j subroutine gives correct outputs.)

using namespace std;

extern "C"
{
    double drc6j_ (double*, double*, double*, double*, double*, 
                double*, double*, double*, int*, int*);

 double w9j_ (int*, int*, int*, int*, int*, 
            int*, int*, int*, int*, int*, int*, double*);

}  
int main()
{
    int N = 1000;
    double* coef9j;
    int L1, L2, L3, L4, L5, L6, L7, L8, L9, Lmin=0, Lmax=6;

    coef9j = new double [ N ];
    L1 = 2, L2 = 2, L3 = 2, L4 = 2, L5 = 2, L6 = 2, L7 = 2, L8 = 2;

    w9j_ ( &L2, &L3, &L4, &L5, &L6, &L7, &L8, &L9, &Lmin, &Lmax, &N, coef9j );

    cout<<" Lmin, Lmax, ierr = "<<Lmin<<","<<Lmax<<","<<ier<<endl;

    for (int i=0;i<Lmax-Lmin+1;i++)
    cout<<"coeff = "<<setprecision(7)<<coef9j[i]<<endl;

    return 0;

}

This program gives values for Lmin and Lmax like E-315 etc. which isn't possible because I should get integer values and that too non zero values for Lmax. In earlier subroutines (like drc6j) the values for Lmin and Lmax were decided by the subroutine, and as I passed by reference, I used those values. Also, it always gives coeff =0 as output regardless of the input Li values.

Please mention in the comments if you would want the drc6j subroutine and the C++ code I use to evaluate it. I'll edit the question accordingly.

The outputs I get for the C++ code are:

AA           2           2           2           2           2           2           2         256           0           6        1000
L9min, L9max, ierr = 0,6,0
coeff = 0
coeff = 0
coeff = 0
coeff = 0
coeff = 0
coeff = 0

And for the Fortran code:

           0           6           7

 main Hi1
  -8.5714285714285719E-003
   1.7347234759768071E-018
   1.6734693877551020E-002
   1.7347234759768071E-018
   1.3877551020408163E-002
   0.0000000000000000     
   0.0000000000000000     
 main Hi3a
 main Hi3b

This is the subroutine w9j code:

subroutine w9j(l1,l2,l3,l4,l5,l6,l7,l8,lmin,lmax,ndim,cof9j)

  implicit none

  !--- in/out variables

  integer, intent(in) :: l1, l2, l3, l4, l5, l6, l7, l8, lmin, lmax, ndim
  real(8), dimension(ndim), intent(out) :: cof9j

  !--- program variables

  integer :: l1min, l1max, l2min, l2max, l3min, l3max
  integer :: i, j, k, jmin, jmax, ier, dim1, dim2, dim3

  real(8) :: y
  real(8), allocatable, dimension(:) :: sixj1, sixj2, sixj3
  real(8) :: rl1, rl2, rl3, rl4, rl5, rl6, rl7, rl8, ri
  real(8) :: rl1min, rl1max, rl2min, rl2max, rl3min, rl3max

  !--- --- ---

  print*, 'AA', l1,l2,l3,l4,l5,l6,l7,l8,lmin,lmax,ndim

  cof9j = 0.0

  do k=1,ndim
    i=k+lmin-1

!    print*, 'i', i
!    print*, ""

    if((l2+i-l1 < 0)  .or. (l2-i+l1 < 0)  .or. (-l2+i+l1 < 0)  .or. &
     & (l2+l5-l8 < 0) .or. (l2-l5+l8 < 0) .or. (-l2+l5+l8 < 0) .or. &
     & (l4+l3-l5 < 0) .or. (l4-l3+l5 < 0) .or. (-l4+l3+l5 < 0) .or. &
     & (l4+l1-l7 < 0) .or. (l4-l1+l7 < 0) .or. (-l4+l1+l7 < 0) .or. &
     & (l6+i-l3 < 0)  .or. (l6-i+l3 < 0)  .or. (-l6+i+l3 < 0)  .or. &
     & (l6+l7-l8 < 0) .or. (l6-l7+l8 < 0) .or. (-l6+l7+l8 < 0)) then
!      print*, "violation of triangularity"
      cof9j(k) = 0.0
      cycle
    end if

    l1min=max(abs(i-l8),abs(l5-l1))
    l1max=min(i+l8,l5+l1)
    dim1 = l1max -l1min + 1
!    print*, 'l1', l1min,l1max,dim1
!    print*, ""

    l2min=max(abs(l3-l7),abs(l1-l5))
    l2max=min(l3+l7,l1+l5)
    dim2 = l2max -l2min + 1
!    print*, 'l2', l2min,l2max,dim2
!    print*, ""

    l3min=max(abs(i-l8),abs(l7-l3))
    l3max=min(i+l8,l7+l3)
    dim3 = l3max -l3min + 1
!    print*, 'l3', l3min,l3max,dim3
!    print*, ""

    if (l1max < l1min .or. l2max < l2min .or. l3max < l3min) then
      cof9j(k) = 0.0
      cycle
    end if

    allocate(sixj1(dim1))
    allocate(sixj2(dim2))
    allocate(sixj3(dim3))

    sixj1 = 0.0
    sixj2 = 0.0
    sixj3 = 0.0

    ri  = real(i)  ! L1
    rl1 = real(l1) ! L2
    rl2 = real(l2) ! L3
    rl3 = real(l3) ! L4
    rl4 = real(l4) ! L5
    rl5 = real(l5) ! L6
    rl6 = real(l6) ! L7
    rl7 = real(l7) ! L8
    rl8 = real(l8) ! L9
    rl1min = 0.0
    rl1max = 0.0
    rl2min = 0.0
    rl2max = 0.0
    rl2min = 0.0
    rl3max = 0.0

!      SUBROUTINE DRC6J (L2, L3, L4, L5, L6, L1MIN, L1MAX, SIXCOF, NDIM, IER)
!          L1MAX=MIN(L2+L3,L5+L6) and L1MIN=MAX(ABS(L2-L3),ABS(L5-L6)).

!    print*, "Hi1"

!    print*, 'BB',ri,rl1,rl2,rl3,rl4,rl5,rl6,rl7,rl8,dim1,dim2,dim3

!              L2  L3  L4  L5  L6
    call DRC6J(ri,rl8,rl2,rl5,rl1,rl1min,rl1max,sixj1,dim1,ier)
!    print*, "6j - 1"
!              L2  L3  L4  L5  L6
    call DRC6J(rl3,rl7,rl4,rl1,rl5,rl2min,rl2max,sixj2,dim2,ier)
!    print*, "6j - 2"
!              L2  L3  L4  L5  L6
    call DRC6J(ri,rl8,rl6,rl7,rl3,rl3min,rl3max,sixj3,dim3,ier)
!    print*, "6j - 3"

!    print*, "Hi2"

    if (l1min /= int(rl1min) .or. l1max /= int(rl1max) .or. &
    &   l2min /= int(rl2min) .or. l2max /= int(rl2max) .or. &
    &   l3min /= int(rl3min) .or. l3max /= int(rl3max)) then
      print*, "Declared dimensions of 'sixj' arrays are wrong"
    end if

    jmin = max(l1min,l2min,l3min)
    jmax = min(l1max,l2max,l3max)
    print*, 'Jmin Jmax', jmin, jmax


    y = 0.0

    if(jmin <= jmax) then

      do j=jmin,jmax
        y = y + (2.0*j+1.0)*sixj1(j-l1min+1)*sixj2(j-l2min+1)*sixj3(j-l3min+1)
      end do

    end if

    cof9j(k) = y

    deallocate(sixj1,sixj2,sixj3)

  end do
  print*, 'end values', l1,l2,l3,l4,l5,l6,l7,l8,lmin,lmax,ndim


  return

end subroutine w9j

The subroutine for drc6j.f and dependencies can be found on netlib.org

Where is the code for the subroutines? Can you show the exact output you get? — Vladimir F Героям слава
– Vladimir F Героям слава, Commented Jun 13, 2016 at 7:43
I don't know if this is important but... in the FORTRAN program you have LL1 LL2... that seem declared as integer, while in th C++ part you declare all parameters as double. This could lead to some problems on the stack, garbling the values of your variables. — Giacomo Degli Esposti
– Giacomo Degli Esposti, Commented Jun 13, 2016 at 7:49
@GiacomoDegliEsposti they seem to have been converted into reals, right? the drc6j subroutine can take only double values (even though here we only compute integers) as you see where it's declared in the extern c function. — pyroscepter
– pyroscepter, Commented Jun 13, 2016 at 7:51
@pyroscepter. You changed to int both the parameters of w9j_ and the variables L1..L9 right? Maybe you have to change also the return type from double to void? — Giacomo Degli Esposti
– Giacomo Degli Esposti, Commented Jun 13, 2016 at 9:05

d_1999 · Accepted Answer · 2016-06-14 11:17:39Z

2

At no point in w9j do you set lmin and lmax and at no point in your c++ code do you set these variables, hence when you print them in your c++ code you are just printing whatever value was in memory at the location of these variables -- i.e. it is undefined. In your original fortran code you do set these variables before passing them to w9j which is why that seems to work.

Note you've marked Lmin and Lmax as intent(in) in w9j but then you're trying to print them out in the c++ code as if they were outputs of w9j. Arguments marked as intent(in) should not be modified within that routine.

Update following comments

The final fix was to note that the c++ call to w9j doesn't pass l1 but you do pass l9 which isn't set to anything, leading to the use of an undefined value in the main w9j routine.

edited Jun 14, 2016 at 11:17

answered Jun 14, 2016 at 8:15

d_1999

8546 silver badges17 bronze badges

Sign up to request clarification or add additional context in comments.

7 Comments

pyroscepter Over a year ago

I modified the code later on and set Lmin Lmax as 0,6 respectively. Still coeff is 0.

pyroscepter Over a year ago

So how do I modify the code to incorporate what you suggested so it doesn't give me 0 as the answer all the time? Why aren't the values being copied to the array coeff, I still have no idea.

d_1999 Over a year ago

@pyroscepter Can you update your question to indicate what code you're actually using now as there have been several suggested changes in the comments. You say you set lmin etc., when you've done this and then print them in your c++ code after the w9j call do you still get unusual numbers?

d_1999 Over a year ago

@pyroscepter Note your c++ call to w9j doesn't pass l1 but you do pass l9 which isn't set to anything.

pyroscepter Over a year ago

That isn't a problem, I had mistaken the code to given an array for l9 values but later I clarified that the array is for l1 values. So my variables here would be shifted by 1, and l9 would be l1. It doesn't matter though since they are all just variables, right?

|

Collectives™ on Stack Overflow

Wrong values with Fortran subroutine called from C++

1 Answer 1

7 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

7 Comments

Your Answer

Sign up or log in

Post as a guest

Related