How to use optional attribute for assumed shape arrays

Question

I encountered the following situation while using the Fortran optional attribute for an assumed shape array and trying to figure out the best (in terms of performance) solution. I would be very glad if someone could give me a good hint. Please note, that I am interested in every performance gain I can get, since my arrays are large and the number of loops and their cycles are even larger.

I have the situation that a calculation is done either using the optional argument or in case it is not present it uses another array at its place.

  subroutine compute(tmat,av,av2)
  implicit none
  complex,intent(out)           :: tmat(:)   
  complex,intent(in)             :: av(:)      
  complex,intent(in),optional    :: av2(:)     

  if(present(av2)) then
      tmat = av *av2
  else
      tmat = av *av
  end if

  end subroutine compute_transition_matrices_bosonic

In this simple example, the above solution would be fine. Nevertheless, my true code is much more complicated. This means that the inner block of the routine may be very large (hundreds of lines) and most important contain many nested loops. One could imagine it like the following:

  if(present(av2)) then
      tmat = "function"(av,av2)
  else
      tmat = "function"(av,av)
  end if

where "function" stands for a lot of operations and loops (therefore the ""). The point is that these are the same operations for both cases of the if statement such that I need to write the code twice. The other solution would be to check if av2 is present at the position of usage in the code which would create some overhead since this check would be done very (very) often.

I was wondering if there is a more clever solution to such a problem. First, one might think of using a temp variable

  complex,                :: phi_tmp(:)

  if(present(av2)) then
      phi_tmp = av2
  else
      phi_tmp = av
  end if

  tmat = "function"(av,phi_tmp)

which is done often when using optional arguments. But this would COPY the data, which is in my case a really huge array. Therefore, I was thinking of using a pointer

  complex,intent(in),target             :: av(:)      
  complex,intent(in),optional,target    :: av2(:)   
  complex,pointer                :: phi_tmp(:)

  if(present(av2)) then
      phi_tmp => av2
  else
      phi_tmp => av
  end if

  tmat = "function"(av,phi_tmp)

But this requires the TARGET attribute for av and av2. Here I am not sure whether this would cause performance drop since the compiler can not assume any more that av and av2 have NO aliases in its optimization procedure, even though here both have the INTENT(IN) attribute such that no aliasing problem can occur. Furthermore, what does it mean if the argument which is in the input when calling the routine does not have the TARGET attribute? (This compiles and works!)

Does anyone has some experience with these issues? Is there a "standard" solution? (I couldn't find one) What ultimately is the best solution?

francescalus francescalus · Accepted Answer · 2015-09-15T10:22:41

One approach would be to use distinct subroutines for the present and not-present cases.

subroutine compute1(tmat,av)
  implicit none
  complex,intent(out)           :: tmat(:)   
  complex,intent(in)            :: av(:)      

  call compute2(tmat, av, av)  ! Ensure there is an explicit interface for compute2 here
end subroutine

subroutine compute2(tmat,av,av2)
  implicit none
  complex,intent(out)           :: tmat(:)   
  complex,intent(in)            :: av(:)      
  complex,intent(in)            :: av2(:)     

  tmat = "function"(av,av2)
end subroutine

where neither has an optional argument.

This can be added to with making compute a generic for both of these.

Naturally, this won't work in all cases, but here with the intent(in) attribute on av and av2 in compute2 we don't have to worry about aliasing.

How to use optional attribute for assumed shape arrays

2 Answers