Assumed-Size Dummy Argument Not Working in Module Procedure

Help
1

I am a bit perplexed by this one, it seems to me that assumed-size dummy arguments are actually treated as rank 1 if the subprogram is in a interface block such as:

randrng subroutine with 2 different interfaces in a module 
module min_random
use, intrinsic :: iso_fortran_env, only: sp => real32
implicit none
private

    interface  randrng
        module procedure randrng_sp
    end interface randrng
    public :: randrng

    public :: randrng_sp

    contains

        subroutine randrng_sp(val, n, val_min, val_max)
        implicit none
            real(sp), intent(out) :: val(*)
            integer, intent(in) :: n
            real(sp), intent(in) :: val_min, val_max
            call random_number(val(1:n))
            val(1:n) = val(1:n)*(val_max - val_min) + val_min
        end subroutine randrng_sp

end module min_random

When I try to write code that calls this function, if I call it as call randrng, then my first argument x seemingly must be an array of rank 1. Otherwise, gfortran, ifort, and ifx all produce errors such as:

ifort compiler error message 
app/main.f90(12): error #6285: There is no matching specific subroutine for this generic subroutine call.   [RANDRNG]
    call randrng(x, size(x), x_min, x_max)
---------^
app/main.f90(16): catastrophic error: Too many errors, exiting
compilation aborted for app/main.f90 (code 1)

The above error is from ifort, but gfortran and ifx say similar.

However, if I call the routine as call randrng_sp then I can supply x of any shape I want, and the first n elements will be filled as originally intended. Is this the way modules and interfaces are meant to work, or why does this happen?

EDIT: For completeness

This program will not compile with gfortran 11.3.0, ifort 2021.9.0, or ifx 2023.1.0 
program main
use min_random
implicit none
    integer, parameter :: n=2
    real :: x(n,n), x_min=0.0, x_max=100.0
    call randrng(x, size(x), x_min, x_max)
    write(*,*) 'x: ',x
end program main
This program will compile and run as expected. 
program main
use min_random
implicit none
    integer, parameter :: n=2
    real :: x(n,n), x_min=0.0, x_max=100.0
    call randrng_sp(x, size(x), x_min, x_max)
    write(*,*) 'x: ',x
end program main
1 Like
9

The direct call works because the storage association between the 2D actual argument and the 1D dummy argument is allowed. The interface access is not allowed because the compiler is matching the type, kind, and rank of the actual and dummy arguments, and the 2D actual argument does not match with the 1D dummy argument.

So to get this to work, you somehow need to get a 1D array aliased to the 2D array, and use that 1D array as the actual argument. There are several ways this can be done, with EQUIVALENCE or with pointers, and maybe some others too. Here is an example that works with pointers.

program main
use min_random
    implicit none
    integer, parameter :: n=2
    real, target :: x(n,n)
    real ::  x_min=0.0, x_max=100.0
    real, pointer :: xp(:)
    call set_xp( x, size(x) )
    call randrng(xp, size(xp), x_min, x_max)
    write(*,*) 'x: ',x
contains
    subroutine set_xp( b, n )
       integer, intent(in) :: n
       real, intent(in), target :: b(n)
       xp => b
       return
    end subroutine set_xp
end program main

There are some situations where a direct pointer assignment like xp=>x will work, but this isn’t one of them, so the assignment is done within the internal set_xp() routine. I think this is standard conforming, so it should be portable.

10

@RonShepard set_xp() should be a function (or have an output argument for the pointer).

11

xp is accessed through host association. Yes, it could also be passed as an argument. I think either approach is standard conforming. A subroutine interface is fine. I generally avoid functions that return pointers because there are so many ways to shoot yourself in the foot, so my first attempt was with a subroutine.

12

(EDITED)

Another possibility is to write a wrapper routine for any array rank, using the assumed rank interface:

module min_random
use, intrinsic :: iso_fortran_env, only: sp => real32
implicit none
private

    interface  randrng
        module procedure randrng_sp_nd
    end interface randrng
    public :: randrng

    contains

        subroutine randrng_sp_nd(val, val_min, val_max)
        implicit none
            real(sp), intent(out) :: val(..)
            real(sp), intent(in) :: val_min, val_max
            select rank(val)
                rank(1) ; call randrng_sp(val,size(val),val_min, val_max)
                rank(2) ; call randrng_sp(val,size(val),val_min, val_max)
                rank default; error stop "not yet written"
            end select
      end subroutine randrng_sp_nd

       subroutine randrng_sp(val, n, val_min, val_max)
        implicit none
            real(sp), intent(out) :: val(*)
            integer, intent(in) :: n
            real(sp), intent(in) :: val_min, val_max
            call random_number(val(1:n))
            val(1:n) = val(1:n)*(val_max - val_min) + val_min
        end subroutine randrng_sp

end module min_random
13

I actually never realized until now that you could query rank for an assumed size array. I had always thought of it as just direct access to memory, in chunks of <variable type kind size> bytes. I see in production code all the time stuff such as: routine is defined only for real of default kind but routine gets called for reals or integers, only the first value of an array is passed to routine with assumed size dummy arguments, etc. All this works because the routines are just free in their own files, and nothing is in modules (which I believe automatically generate interfaces for them).

Or at least, I thought that putting code in modules automatically generated an interface for them. I still don’t understand why the direct call (still with routine defined in a module) is allowed while the interface one is not. Your post says “the storage association is allowed,” but I do not understand what that means. Why does the direct call not complain about mismatched rank given assumed size dummy argument? It would if the dummy argument was assumed shape.

I don’t love any solution that is ultimately “write a wrapper routine and call it that way,” but alas it seems that is maybe the only way at the moment. I absolutely do not want to use assumed rank because it does not work correctly in gfortran at the moment.

14

Note the difference between v(..) and v(*) in @PierU’s example. The former is assumed-rank, the latter assumed-size.

15

I see that you are quoting text from multiple different posts. Thank you for your reply, and yes I understand the difference between assumed rank and assumed size. I refuse to use assumed rank for the time being due to lack of proper support in gfortran.

I have already marked the quoted post from @RonShepard as solution, but for my own understanding - I do not know what “storage association” means. As I posted above:

1 Like
16

Sorry, I missed that paragraph. Storage association is Fortran standard jargon for associating a variable with a specific storage location in memory.

The direct call

real(sp) :: x(n,n)
call randrng_sp(x,size(x),x_min,x_max)

behaves as if only the address of the first element in x is passed.

We could compare this to array decay in C, i.e.

void randrng_sp(float val[], int *n, float *val_min, float *val_max);

is equivalent to

void randrng_sp(float *val, int *n, float *val_min, float *val_max);

Hence, why you can use assumed-size arrays in C bindings, because the concept maps closely.

17

Thank you very much that was a clear and concise explanation. I see this type of assumption made all the time in older code. Often calls to subroutines taking assumed size arrays are made and the caller only provides the element from which they wish the callee to begin.

It is a shame that these same rules are not applied to code grouped within interface blocks in modules. It is a very nice functionality.

18

I missed it…

Well, you can’t, actually. Inquiring the rank is possible only on an assumed rank array.

It does not complain because associating an actual argument to an assumed size and or an explicit shape dummy argument of a different rank is explicitely allowed by the standard. But as you have pointed out in this discussion, it’s not possible if the routine is called throught a generic interface.

I compiled and run my example with gfortran. I don’t know what is the problem with assumed rank in gfortran, but I don’t think this kind of simple code may break.

19

There are often some reasons why some given feature is not available in some context.

20

Actually you can use iso_c_binding and write:

subroutine randrng_sp_nd(val, val_min, val_max)
use iso_c_binding
implicit none
    real(sp), intent(out),contiguous,target :: val(..)
    real(sp), intent(in) :: val_min, val_max
    real(sp), pointer :: vp(:)
    
    call c_f_pointer(c_loc(val), vp, [size(val)])
    
    call randrng_sp(vp,size(val),val_min, val_max)
end subroutine randrng_sp_nd

And that will work for every rank.
I would have liked to write:

vp(1:size(val)) => val

But it seems it isn’t possible…

21

Note the Fortran standard committee, or at least a member will be inclined to remind you of a line from Shakespeare! If not that, then this proverb!

Your issue is the generic interface. You make the specific procedure randrng_sp accessible as PUBLIC anyway, so have your caller consume that rather than the generic name. With assumed-size received argument in that specific procedure, the generic interface is meaningless as it can be ambiguous given the overall semantics laid out in the language.

Assumed-size received arguments can be convenient but they have aspects to them that can be like chainsaws without bladeguards. Generic interfaces are like newer “tools” that strive for some safety. Avoiding a mix-and-match of the two is something you can consider until Generics facility in the language with better support toward RANK genericity come about a decade or two from now.

22

Excuse me, but I disagree. In modern Fortran, you can write code that does what it look like it does, in the spirit of FORmula TRANslation. If you want to pass x(3:7) you do that instead of passing arguments x(3) and n=5. If your arrays are 3-D, you declare them as 3-D and not 1-D and pass them as such. If you want to pass a 3-D array section as 1-D, you can use reshape or the feature that [x] is a 1-D array regardless of the rank of x. I don’t think new code should be written that relies on storage association.

23

To the extent OP is open at all to the idea of using an additional subprogram a la a helper procedure, OP can also consider a shorter option with the internal procedures introduced in the language starting Fortran 90:

program main
    implicit none
    integer, parameter :: n=2
    real, target :: x(n,n)
    real ::  x_min=0.0, x_max=100.0
    call helper( x, size(x) )
    write(*,*) 'x: ',x
contains
    subroutine helper( a, n )
       use min_random
       real, intent(inout) :: a(*)
       integer, intent(in) :: n
       call randrng(a, n, x_min, x_max)
       return
    end subroutine
end program main
24

Well, in turn I partly disagree here. I indeed see no reason to pass x(3) if your intent is to pass x(3:7), but storage association is still useful when it comes to pass an array with a rank different from the rank declared in the routine. The reshape() solution is not always a good one: not only it can result in a temporary array, but since it is an expression it cannot be used if the dummy argument is intent(out).

The assumed rank feature is a possible approach, but as shown in my example at the end the deepest call still relies on storage association.

The recent versions of the standard have introduced the possibility to point to a 1D array with a nD pointer, and I think the opposite, i.e. pointing to a nD array with a 1D pointer, would be useful. Of course this probably assumes that the target is contiguous.

25

It’s possible, and the target doesn’t have to be contiguous:

! test_remap.f90
program remap
implicit none
real, target :: a(2,2)
real, pointer :: b(:)
a = reshape([1,2,3,4],[2,2])
b => a(2,:) ! second row [2,4]
print *, b 
print *, is_contiguous(b)
end program

~/fortran$ gfortran -Wall test_remap.f90 
~/fortran$ ./a.out
   2.00000000       4.00000000    
 F
26

When I said the storage association is allowed, I meant that the fortran standard defines what a program does that uses storage sequence association. It is allowed by the standard, its behavior is defined, and it therefore should be portable. This does not depend on byte order (i.e. big/little-endian) conventions, it is more portable than that. In fact it would work correctly even on a word-addressable machine, with no bytes anywhere.

As you say, this storage association would not be allowed if the dummy argument were assumed shape. So if you want to use storage sequence association, then don’t do that. In addition to assumed-size arrays as the dummy argument, explicit shape arrays could also be used.

A downside of this approach is that it places all of the bounds checking responsibility on the programmer. The compiler cannot help if you pass in, for example, a four element array and then reference six elements within the subroutine. Another problem, at least performance-wise, is that if you pass in a noncontiguous actual argument, say an array slice of some kind, or an array of derived type components, then the compiler will be required to make a temporary copy in order to associate the actual and dummy arguments. So that is allowed, and supported, by the language, but the copy-in/copy-out may result in a big performance hit, and it may not be obvious why your code is suddenly running 10x slower without a lot of debugging effort.

In my earlier post, I said that you needed to somehow alias a 1D and 2D array. Another way of doing that was posted by egio.

This has the advantage that it does not require the user to write something like subroutine set_xp( ), but it is doing exactly the same thing inside, it is making the 1D array an alias to the 2D array. So for purposes of demonstrating the technique, a set_xp() subroutine might be best, but in practice it might be easier to just use the intrinsic subprogram. The only problem might be that your fortran compiler does not support the intrinsic module iso_c_binding, in which case you would need to fall back and write your own version. Both approaches are relying on the same underlying storage sequence association.

1 Like
27

I meant pointing to a whole nD array with a 1D pointer. c_f_pointer(c_loc(a),b,[size(a)]) does the trick, but it artificially relies on C binding.