In F90 we can have functions that return array results. If the input arguments are assumed-shape, the result must be automatic. For example, we can reformulate the commutator as an array-valued function:
function commut(A,B) implicit none real, dimension(:,:), intent(in) :: A,B real, dimension(size(a,1),size(a,2)) :: COMMUT if( (size(A,1)/=size(A,2)) .or. & (size(A,2)/=size(B,1)) .or. & (size(B,1)/=size(B,2)) ) then print*, 'Error: input matrices not conformal' return end if commut = matmul(A,B)-matmul(B,A) end function commutHere we check the compatibility of the input arguments.
Note that declaring the result (COMMUT) as an assumed-shape array does not work, since it is not associated with any actual argument; we have to declare COMMUT explicitly with proper dimensions that depend on the size of dummy arguments - in other words, COMMUT is declared as an automatic array.