Opaque Data

[Dishaw__Jim]

I have the following problem that I am trying to solve. I have a subroutine that implements the spatial geometry independentportion of the algorithm. One of the arguments that is passed to that subroutine is a geometry specific subroutine, which implements the geometry specific part of the algorithm.

Driver Program ----> Solver(geometry_solver) ----> geometry_solver()

I had everything working when I made the result arrays sized for a 3D problem (2D and 1D problems had the extraneous dimensions sized to 1). While this approach does work, I find it a bit crude in the sense that I am making a rank 3 array when only a rank 1 or 2 is needed.One other approach that I can think of is to make a phi_1d, phi_2d, and phi_3d array and use the appropriate one. This solution just shifts the kludge elsewhere.

An alternative approach is to make a rank 1 arrayand eithermanually index the array or useRESHAPE. From what I understand RESHAPE will copy the data. I would prefer not to manually index the array because that obscures what the code is doing.

In C, I could have a pointer to an opaque piece of data and the pointer gets casted when it needs to be used. I looked at the integer pointer, however, that is not standard Fortran.

My preference is to stay with standard Fortran 95 or 2003. Any recommendations?

[Steven_L_Intel1]

One of my favorite bits from the C Interoperability features of F2003 (supported in 10.1) is the ability to do the equivalent of casts using the C_PTR type and support for it. Let's say you have an array that is rank 3 that you want to "cast" to a rank 1 array. You can do something like this:

real, pointer, dimension(:) :: rank1_array
...
call c_f_pointer (c_loc(rank3_array), rank1_array, [size(rank3_array)])

Now rank1_array is a rank 1 array of the proper dimensions and no copying. Here's how it works. c_loc returns a "c pointer" to its argument, rank3_array. That value is a derived type C_PTR. The subroutine c_f_pointer converts a C_PTR to a Fortran pointer (second argument), assigning it the shape of the third argument, which we have specified as a rank 1 array containing one element, the number of elements in rank3_array.

If you wanted to extend this to a rank 2 pointer, you'd declare:

real, pointer, dimension(:,:) :: rank2_array

and then make the third argument to c_f_pointer be an array of two elements being the number of elements in each dimension. (I leave that as an exercise to the reader.)

An alternative, which might or might not be better for you, is to take advantage of Fortran "sequence association", where you can pass the first element of the array and it can be passed to an array argument of any rank.

[Dishaw__Jim]

MADsblionel:

An alternative, which might or might not be better for you, is to take advantage of Fortran "sequence association", where you can pass the first element of the array and it can be passed to an array argument of any rank.

Thanks for the tip. I always forget about the "sequence association" because I always feel like I am breaking the rules. Unfortunately, I do not think it will work in my case because I get the following error: "If dummy arg is allocatable, actual arg must be a whole array and not a section." Below is an example of what I am trying to do:

! 1D example
SUBROUTINE foo
  REAL, ALLOCATABLE :: data(:)

  CALL setup(1, data)
END SUBROUTINE foo

SUBROUTINE setup(geometry, data)
  REAL, ALLOCATABLE, INTENT(OUT) :: data(:)

  SELECT CASE(geometry)
  CASE (1)
    CALL SETUP_1D(data)
  CASE (2)
    CALL SETUP_2D(data)
  CASE (3)
    CALL SETUP_3D(data)
  END SELECT
END SUBROUTINE

SUBROUTINE seutp_1D(data)
  REAL, ALLOCATABLE, INTENT(OUT) :: data(:)

  ! 1D setup here
END SUBROUTINE setup_D

SUBROUTINE seutp_2D(data)
  REAL, ALLOCATABLE, INTENT(OUT) :: data(:,:)

  ! 2D setup here
END SUBROUTINE setup_2D

SUBROUTINE seutp_3D(data)
  REAL, ALLOCATABLE, INTENT(OUT) :: data(:,:,:)

  ! 3D setup here
END SUBROUTINE setup_3D

I guess the C interoperability feature might be the way to go. I haven't tried using POINTER, but I don't think it will work either.

[Dishaw__Jim]

Incidentally, part of the reason I have not pursued using POINTER is that you had mentioned the compiler will forego some optimizations because it does not know if the memory is contiguous.

[Steven_L_Intel1]

Without seeing more of your code I can't tell for sure, but I might guess that you don't need or want ALLOCATABLE in those dummy argument declarations, unless you're planning to ALLOCATE or DEALLOCATE them in the routines.

I believe I said that I doubted you would be able to tell the difference in a real application. My advice is to use the language features the way they were intended to be used and to not avoid a feature just because you fear a loss of optimization - especially if the alternative greatly complicates the code. The compiler does well with POINTERs in most cases.