Fortran Function: best practice when output is related to implicit size input

Question

For example, I want a function to calculate the commutator of two matrices. The way I implemented it is like this:

   function commutator_dd(MA,MB)
 > result(MC)
   implicit none
   real*8,intent(in)      :: MA(:,:)
   real*8,intent(in)      :: MB(:,:)
   real*8,allocatable     :: MC(:,:)
   real*8,allocatable     :: MP(:,:),MN(:,:)
   integer                :: nn

   nn=size(MA,1)
   allocate(MC(nn,nn),MP(nn,nn),MN(nn,nn))
   MP=MATMUL(MA,MB)
   MN=MATMUL(MB,MA)
   MC=MP-MN
   return;end function

A friend of mine pointed out that perhaps allocating arrays and not deallocating them was not a good idea. I would tend to agree with that, but as it is I cannot deallocate MC before returning to main (I think fortran automatically deallocates them after exiting the procedure, but from what I remember leaving things like this was not deemed a good practice). He suggested giving nn as an input to the function and declaring all matrices as M(nn,nn) instead of using allocation (or at least doing it with MC), but I kind of dislike the idea of adding unnecessary inputs.

Other possibility would be to declare MC(size(MA,1),size(MA,1)), but I had issues with this sort of syntax before (I can't pinpoint what exactly the problem was because It was always easier to circumvent it by using a different method of passing the variable, but I remember getting complains from the compiler when trying to do it that way).

Which would be the best way to do this? which is the way in which the intrinsic procedures (like matmul) deal with this? (I tried to google for matmul's code but what I found was definitively not as straightforward as I was expecting).

Answer 1

The automatic array (I am not sure it is the proper standard term for a function result) which @HighPerformanceMark suggest are the best design choice here, but you should not worry even about the allocatable arrays.

In modern Fortran any memory leak is simply impossible with anything allocatable being it a scalar or an array, be it a variable or a derived type component, whatever. The memory is always guaranteed to be deallocated automatically, unless it has the save attribute (which you normally use for that purpose).

In your case the allocatable function results are automatically deallocated by the calling code, while the temporary arrays are automatically deallocated when the subroutine finishes.

Also note that in Fortran 2003 (which you need for the allocatable function results anyway) you can omit the allocate statement. The arrays will be allocated automatically during the = assignment.

Summary: You rarely need deallocate() for anything allocatable. You use it more often for pointer variables. Regarding what is best practice - why clutter the code with unnecessary statements, when the correct behaviour is guaranteed by the standard?

Answer 2

Me, I'd just use an automatic array for the return, like this

function commutator_dd(MA,MB) result(MC)
   implicit none
   real*8,intent(in)      :: MA(:,:)
   real*8,intent(in)      :: MB(:,:)
   real*8                 :: MC(size(ma,1),size(ma,2))

   MC=MATMUL(MA,MB)-MATMUL(MB,MA)
   return
end function

and let the compiler take care of creating temporaries if it wants to. Like you I see no need (not since about Fortran 90) for including the dimension(s) of the input arrays in the argument list. Of course this snippet does no error checking and fails ungracefully if the input arrays are not conformable.

If your compiler does complain about the automatic sizing of MC report back; I can't see any reason why it should.