Messing with POINTERS and ALLOCATABLE dimensions

[holysword]

I'm working with this code which basically allocates and deals with arrays of multiple kinds and dimensions, e.g.

SUBROUTINE calc_real_2(u)
   REAL, ALLOCATABLE, INTENT(INOUT) DIMENSION(nx,ny) :: u
.
.
.

SUBROUTINE calc_real_3(u)
   REAL, ALLOCATABLE, INTENT(INOUT), DIMENSION(nx,ny,nz) :: u
.
.
.

SUBROUTINE calc_cmplx_3(u)
   COMPLEX, ALLOCATABLE, INTENT(INOUT), DIMENSION(nx,ny,nz) :: u
.
.
.

They are all provided in an interface which gathers them all into one "calc" function. The user is supposed to pass u and the rest is decided by the types involved. Since I have to deal with reals integers and complex numbers, ranging from 2 to 7 dimensions, I'd have an ridiculously large number of functions doing basically the same thing. You know, one of those situations you wish you had C's templates.

I'm wondering if it is possible to simplify this. I could, of course, have one function for each type looking like

SUBROUTINE calc_real(u)
   REAL, ALLOCATABLE, INTENT(INOUT), DIMENSION(nx*ny*nz) :: u
.
.
.

and so on, treating the indexes lexicographically. If it is 2D, then nz=1 and everything should work fine. The indexes would look dreadfully complicated for 5D or 6D though, and elegance is also something important if we are expecting to maintain a code. I'd like to use EQUIVALENCE, but it doesn't work for pointers or dummy arguments or allocatables or anything dynamically allocated (which cuts its usefulness by half, really). Maybe allocating another variable and using RESHAPE, but that would copy the contents from one to another, right?

Besides, having it declared as u(:) kinda destroys the overloading in this case. If I have w(:,:,:) declared somewhere and I call calc_real(w), it works fine, the compiler seem to force w to be interpreted as a 1D array. However if I use the interface "calc",  call calc(w) says that no matching function was found. I've considered using the (*) but that's usually hurts the performance if I am not mistaken.

I'm sure some of you have stumbled upon this kind of problem once or twice, how did you walk around it?

[Steven_L_Intel1 · Accepted solution]

What you want is "assumed rank", a Fortran 2015 feature we'll support in the 16.0 compiler (beta in progress.) You would declare U as DIMENSION(..) and could pass any rank argument to it. Once you get into the routine, you'd have to use C_LOC and C_F_POINTER to "cast" it to an array of the proper dimension, or to one that is rank 1 and you do the index computation, or whatever.  There is a proposal for a "SELECT RANK" construct to make this easier, but it hasn't yet been accepted into the language.

I'll comment that you need more than 7 dimensions - the standard currently specifies a minimum of 15 and we support up to 31.

[FortranFan]

Steve Lionel (Intel) wrote:

What you want is "assumed rank", a Fortran 2015 feature we'll support in the 16.0 compiler (beta in progress.) You would declare U as DIMENSION(..) and could pass any rank argument to it. Once you get into the routine, you'd have to use C_LOC and C_F_POINTER to "cast" it to an array of the proper dimension, or to one that is rank 1 and you do the index computation, or whatever.  There is a proposal for a "SELECT RANK" construct to make this easier, but it hasn't yet been accepted into the language.

..

Steve,

Good point about the new assumed rank feature.  Have you , by any chance, tried out some code using the new feature in 2015 standard with the 16.0 beta compiler that you can share here?

OP mentioned C++ templates implying generic programming.  Unfortunately Fortran still falls short by some measure in this aspect.  I wish the parameterized derived types (PDT) feature would get enhanced further in some not-too-distant future to support more generic programming features.  Something like the following will be cool if it were to become part of the standard!

module m
   
   implicit none
   
   !.. pseudo code
   type, public :: foo(T) !.. Similar syntax as PDT
      
      generic, type :: T  !.. similar to <integer, kind> syntax in PDT
      
      type(T), dimension(..), allocatable :: u  !.. similar to assumed rank and PDT syntax
      
   contains
   
      procedure, pass(this) :: calc_u
      
   end type foo
   
contains

   subroutine calc_u(this, ..)
      
      class(foo(T)) :: this
      
      
      .. 
      select type ( this%u )
      
         type is (integer)

            ..  !.. may use <select rank> feature here if it gets added

         type is ( real )
            ..
         type is ( complex )
            ..
      end select
      
      .. !. for many common operations, select type construct may not be needed!
      
   end subroutine calc_u
   
end module m

program p
   
   use m, only : foo
   
   type(foo(integer)) :: a  !.. integer type for u
   type(foo(real))    :: b  !.. real type for u
   
   allocate( a%u(0:2), source=[1,2,3] )  !.. here u is 1-D array with 3 elements
    
   call a%calc(..)  !.. some calc on u
   
   allocate(b%u(-10:10,1:3,0:2))  !.. here u has rank 3 with the supplied shapes!
   
   ..
   
end program p

Something you can perhaps add to your notes to bring to a standard committee meeting in the future!  I can see something like above being very useful for many container classes  - lists, collections, stacks, hashsets, queues, etc in addition to the problem posed in the original post here.

P.S.> Welcome back, by the way - hope you'd a great time-off!

[Steven_L_Intel1]

Your application seems appropriate for CLASS(*) - you can do SELECT TYPE on it. Yes, I have tried out pretty much all of the TS29113 features we added, including all of the ISO_Fortran_binding.h procedures.

Here's an example from my webinar of two weeks ago (still waiting for the recording...) It doesn't use dimension(..) though.

use, intrinsic :: iso_c_binding
interface
  function c_alloc (array) bind(C)
  import
  integer(C_INT) :: c_alloc
  real(C_FLOAT), intent(out), allocatable, dimension(:) :: array
  end function c_alloc
end interface

real(C_FLOAT), allocatable, dimension(:) :: my_array
if (c_alloc(my_array) == 0) then
  print *, lbound(my_array), ubound(my_array); print *, my_array
  end if
end

#include "ISO_Fortran_binding.h“
extern int c_alloc (CFI_cdesc_t * descr) {
	int ret, i; float * array;
	CFI_index_t lower = 0, upper = 10;
	ret = CFI_allocate (descr, &lower, &upper, 0); // No elem_len
	if (ret == CFI_SUCCESS) {
		array = descr->base_addr;
		for (i=lower;i<=upper;i++) {array = (float) i;}
	}
	return ret;
}

C:\Projects>icl -c c_sub.c
Intel(R) C++ Intel(R) 64 Compiler for applications running on Intel(R) 64, Versi
on 16.0.0.049 Beta Build 20150501
Copyright (C) 1985-2015 Intel Corporation.  All rights reserved.

c_sub.c
icl: NOTE: The Beta evaluation period for this product ends on 25-sep-2015 UTC.

C:\Projects>ifort fmain.f90 c_sub.obj
Intel(R) Visual Fortran Intel(R) 64 Compiler for applications running on Intel(R
) 64, Version 16.0.0.049 Beta Build 20150501
Copyright (C) 1985-2015 Intel Corporation.  All rights reserved.

ifort: NOTE: The Beta evaluation period for this product ends on 25-sep-2015 UTC
.
Microsoft (R) Incremental Linker Version 12.00.31101.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:fmain.exe
-subsystem:console
fmain.obj
c_sub.obj

C:\Projects>fmain.exe
           0          10
  0.0000000E+00   1.000000       2.000000       3.000000       4.000000
   5.000000       6.000000       7.000000       8.000000       9.000000
   10.00000

[FortranFan]

Steve Lionel (Intel) wrote:

Your application seems appropriate for CLASS(*) - you can do SELECT TYPE on it. ..

Per current standard, yes unlimited polymorphic object is the only option for such needs.  But as you know, that provides little to no compile-time benefits.  In addition, the "consumer" of the class too has to "unbox" the object using <select type> constructs, a significant burden.  Something like above can avoid such issues.

[Steven_L_Intel1]

I have seen suggestions of templates before. It may be that a proposal comes up for the next revision. In your example, I guess the compiler would have to invent variations on calc_u, one for each possible type of T, unless you wanted it done at runtime, in which case why not use class(*)?

[FortranFan]

Steve Lionel (Intel) wrote:

... In your example, I guess the compiler would have to invent variations on calc_u, one for each possible type of T, ..?

Yes, very much so - I'm trying to pass the burden onto the compiler!  But note no more than what compilers of other object-oriented languages do.  You'll notice the example I gave above has the same general concept as a generic container class in Java:

class Container<T> { 
  T m_ob; 

  Container(T o) { 
    m_ob = o; 
  } 

  T getob() { 
    return m_ob; 
  } 

} 

And one can do the same with several other OOP languages too e.g. ,Microsoft seems to have copied shamelessly from Java for its .NET languages.  I'm hoping one can bring similar capabilities to Fortran.

[Steven_L_Intel1]

Ok, but what one does in the other languages is instantiate an object of the class with the desired parameters, typically using "new". You then reference everything through the class object. In a Fortran sense, this is closer to type-bound procedures.

I'm sure there is a proposal floating around somewhere and it is likely to be considered for the next revision (F2020?)

[holysword]

I went out for the weekend and then came back to this whole lot of information. Thank you very much!

I was wondering "where is Steve Lionel, he always replies within 2 or 3h..." so there you have it, he was not around, that explains everything.

Good to know that there is a proposal for SELECT RANK. I could play with those C_LOC and C_F_POINTER now, but I wouldn't dare to rely on it for the production code yet due to portability issues.

 In your example, I guess the compiler would have to invent variations on calc_u, one for each possible type of T

I don't know about FortranFan, but that's what I would like to see, specially when the only thing changing is the KIND of the intrinsic types involved.