1

To put some background context, I need a thread safe random number generator for use in multithreaded Fortran code that needs to be cross compiler, and cross platform compatible. The best way to achieve this is to stick with the language standards. So I wanted to wrap the C++11 random number generators in a function I can call from fortran so I have one random number generator per thread each with it's own state object.

So I created a small C++ class with 3 C function bindings.

#include <random>
#include <iostream>

class random {
    std::mt19937_64 engine;
    std::uniform_real_distribution<double> dist;

    public:
    random(uint64_t seed) : engine(seed), dist(0.0, 1.0) {};
    double get_number() {
        return dist(engine);
    }
};

extern "C" {
    void *random_construct(int seed) {
        return new class random(static_cast<uint64_t> (seed));
    }

    double random_get_number(void *r) {
        return static_cast<class random *> (r)->get_number();
    }

    void random_destroy(void *r) {
        delete static_cast<class random *> (r);
    }
}

A Fortran interface

  MODULE random
  USE, INTRINSIC :: iso_c_binding

  IMPLICIT NONE

  INTERFACE
     TYPE (C_PTR) FUNCTION random_construct(seed)                          &
 &   BIND(C, NAME='random_construct')
     USE, INTRINSIC :: iso_c_binding

     IMPLICIT NONE

     INTEGER (C_INT), VALUE :: seed

     END FUNCTION
  END INTERFACE

  INTERFACE
     REAL (C_DOUBLE) FUNCTION random_get_number(r)                         &
 &   BIND(C, NAME='random_get_number')
     USE, INTRINSIC :: iso_c_binding

     IMPLICIT NONE

     TYPE (C_PTR), VALUE :: r

     END FUNCTION
  END INTERFACE

  INTERFACE
     SUBROUTINE random_destroy(r)                                          &
 &   BIND(C, NAME='random_destroy')
     USE, INTRINSIC :: iso_c_binding

     IMPLICIT NONE

     TYPE (C_PTR), VALUE :: r

     END SUBROUTINE
  END INTERFACE

  END MODULE

And a small program to test this.

   PROGRAM random_test
   USE random

   IMPLICIT NONE

   TYPE (C_PTR) :: rand_object
   INTEGER      :: count

   CALL SYSTEM_CLOCK(count)
   rand_object = random_construct(count)

   WRITE(*,*) random_get_number(rand_object)

   CALL random_destroy(rand_object)

   WRITE(*,*) random_get_number(rand_object)  ! Expected to segfault.

   END PROGRAM

Running this shows my destroy function does not appear to be working correctly since calls after my destroy function are still generating random numbers.

If I change my test program to

   PROGRAM random_test
   USE random

   IMPLICIT NONE

   TYPE (C_PTR), ALLOCATABLE :: rand_object
   INTEGER                   :: count

   CALL SYSTEM_CLOCK(count)
   rand_object = random_construct(count)

   WRITE(*,*) random_get_number(rand_object)

   DEALLOCATE (rand_object)

   WRITE(*,*) random_get_number(rand_object)

   END PROGRAM

Now it's producing the behavior I would have expected and segfaults after the DEALLOCATE. Now my gut reaction says this shouldn't work since I would never try to allocate memory in one language and deallocate it in another. But is there any reason why it shouldn't? The C++ object is a POD type so its memory should be continuous. Then as long as Fortran has the right memory address, it should be able to trivially deallocate it. Is there something I'm missing here?

user1139069
  • 1,505
  • 2
  • 15
  • 27
  • 2
    Accessing a deleted object is undefined behavior. It's not required to crash. – François Andrieux Sep 18 '18 at 17:54
  • Except if I explicitly, null the object after calling delete, 'void random_destroy(void *r) {delete static_cast (r); r = nullptr;}' It still doesn't generates random numbers. – user1139069 Sep 18 '18 at 17:57
  • That won't change anything. `r` is a copy of the pointer and assigning another address to it won't affect the original pointer in any way. – François Andrieux Sep 18 '18 at 17:58
  • That doesn't make any sense. You cannot call a method on a memory address of 0. – user1139069 Sep 18 '18 at 17:59
  • Adding `r = nullptr;` to `random_destroy` will not assign null to `rand_object`. `r` is a copy of `rand_object`, so assigning to `r` doesn't change `rand_object`. – François Andrieux Sep 18 '18 at 18:01
  • Fortran is pass by reference. Why shouldn't it? – user1139069 Sep 18 '18 at 18:05
  • `r` is a reference to the object that `rand_object` refers to yes, in the sense that they both refer to the same object. But `r` is not a reference to the reference `rand_object` itself. Edit : I'm not an expert in Fortran, but that's my understanding, and it's consistent with what you describe. If `rand_object` was nulled, then `random_get_number` should in all likelihood crash. – François Andrieux Sep 18 '18 at 18:15
  • `TYPE (C_PTR), ALLOCATABLE :: rand_object` declares an allocatable scalar C pointer. This is nothing to do with an array object/container. – francescalus Sep 18 '18 at 18:24
  • 1
    @user1139069 I think it is because you are attaching the VALUE attribute in your interface for the external routines. If you remove VALUE and change the definition of C++ routines such that they receive void **r (and setting *r = nullptr in the destroy routine), it may work... (I've tried it on my mac, and it seems to be working.) – roygvib Sep 18 '18 at 18:25
  • @user1139069 [But I remember there is a good thread-safe MT PRNG in Fortran somewhere on the net, so it may be easier to use it (if available...)] – roygvib Sep 18 '18 at 18:27
  • @roygvib That's what I was missing. I forgot that the C side arguments would show up as a pointer to the argument. `TYPE(C_PTR)` is the equivalent of `struct {void *}` explaining why my fortran and C addresses kept coming up different. Unfortunately, someone can no longer add a correct answer since it was erroneously marked as a duplicate. – user1139069 Sep 18 '18 at 19:57
  • I won't vote yet to reopen this question. Despite the title, there's nothing in the question to suggest that you are deallocating memory in Fortran that was allocated in C++. Regarding the `value` attribute: the C++ routine which does the destruction clearly can't update the Fortran argument which was passed by value. This aspect is covered in the new question I added to the duplicate list. – francescalus Sep 18 '18 at 22:02
  • That's the point. Did you miss the "Is there something I'm missing here? " – user1139069 Sep 19 '18 at 18:02
  • Assuming the last comment was addressed to me: you have missed something, but that something is covered in the other question. Please be assured that having a question marked as a duplicate isn't a statement that the question is useless in some way. It's still a valid question you have. But if you think the other question _doesn't_ apply, then you can say so. – francescalus Sep 19 '18 at 19:31

0 Answers0