1 Topic by Boris_MV (edited by Boris_MV 2015-01-13 00:58:02)

Topic: Problem declaring function which outputs vector

Hi there,

I have following issue with declaring a very simple function (CrossProduct) which outputs a vector as a result (which is cross product of input vectors a and b), this error is reported for the declaration line:  "EXTERNAL attribute conflicts with DIMENSION attribute."
So, it tells me I can't do declaration (see the code below):
REAL(kind=8), DIMENSION(3), EXTERNAL :: CrossProduct

Here is the excerpt from the code:

SUBROUTINE Rut1(m,n, v)
    IMPLICIT NONE
    REAL(kind=8), DIMENSION(3), EXTERNAL :: CrossProduct
        ...
    ...
END SUBROUTINE Rut1

!--------------------------------------

FUNCTION CrossProduct(a,b)
    IMPLICIT NONE
    REAL(kind=8)   :: a(1:3), b(1:3)
    CrossProduct(1) = a(2) * b(3) - a(3) * b(2)
    CrossProduct(2) = a(3) * b(1) - a(1) * b(3)
    CrossProduct(3) = a(1) * b(2) - a(2) * b(1)
END FUNCTION CrossProduct

THANK YOU VERY MUCH IN ADVANCE.

Re: Problem declaring function which outputs vector

Fortran doesn't allow you to use functions which return arrays unless there is an explicit interface. EXTERNAL doesn't create an explicit interface so the use of this attribute and the DIMENSION attribute together is forbidden.

Either put the CrossProduct in a module or use an INTERFACE block.

For example, the following works OK. I added the use of INTENT(in) for completeness but this isn't strictly necessary.

SUBROUTINE Rut1
   IMPLICIT NONE
   
   real(kind=8), dimension(3) :: a, b, c
   
   interface
      function CrossProduct(a,b)
         real(kind=8), intent(in), dimension(3) :: a, b
         real(kind=8), dimension(3) :: CrossProduct
      end function
   end interface
   
   a = (/ 1, 0, 0/)
   b = (/ 0, 1, 0/)
   c = CrossProduct(a, b)
   print *, c

END SUBROUTINE Rut1

!--------------------------------------

FUNCTION CrossProduct(a,b)
    IMPLICIT NONE
    REAL(kind=8), intent(in), dimension(3)   :: a, b
    real(kind=8), dimension(3) :: CrossProduct
    CrossProduct(1) = a(2) * b(3) - a(3) * b(2)
    CrossProduct(2) = a(3) * b(1) - a(1) * b(3)
    CrossProduct(3) = a(1) * b(2) - a(2) * b(1)
END FUNCTION CrossProduct

! Test program
program anon
   call Rut1
end
--
David

3 Reply by Boris_MV (edited by Boris_MV 2015-01-14 01:25:38)

Re: Problem declaring function which outputs vector

Wonderful explanation. Thanks a lot David!

Re: Problem declaring function which outputs vector

Interface blocks didn't exist whilst some of us were learning about spaghetti, with Fortran IV, Fortran 66 and Fortran 77.

I have managed to lose most of my worst bad habits since then, but not all.
I've stopped dimensioning arrays as size (1) then deliberately exceeding array bounds.  I greatly prefer "REAL*8" to "REAL (KIND=8)" and I am still finding it unnatural to use double-colons. 

Instead of 

real(kind=8), dimension(3) :: a, b, c

I still write 

REAL*8 a(3), b(3), c(3)

I've almost stopped using INCLUDE statements and COMMON blocks.  I also do now use MODULES, but mainly as a substitute for INCLUDE statements and argument lists, so I am sure I'm using them very inefficiently.  And I have not yet grasped what exactly an interface block does.  Is it analogous to an argument list in a call statement?  If yes, then does it have to be repeated in each program unit that uses the module?

I do intend to take the trouble to read up on this, but any brief hints from others with experience of migrating to modern Fortran would be appreciated.

5 Reply by davidb (edited by davidb 2015-01-14 22:13:20)

Re: Problem declaring function which outputs vector

Some subroutines and functions must have an explicit interface for them to be called with certain arguments or to return certain function results. Examples are array functions (like above), recursive functions, assumed shape array arguments.

An interface block is just a way for the programmer to tell the compiler how a subroutine or function is intended to be called and how the actual arguments are associated to the dummy arguments. For example, for assumed shape arrays information about the size of the array and stride length are also passed (secretly).

The interface block also allows the compiler to check that the calls made actually agree with the interface (argument list, any intents, and function type). An error will be generated if there is a mismatch in any of the arguments with the interface.

However, using interface blocks in this way isn't really good Fortran because:

(1) You have to manually copy/paste the interface block to each place where you want to make the call, and
(2) The compiler isn't actually checking that the subroutine being called has the declared interface.

A much better way to tell the compiler about the interface is to put the subroutine or function in a module. This automatically gets around the above limitations.

You should use an interface where the function or subroutine cannot be placed in a module. Perhaps it is old F66 or F77 code you don't want to change, some functions in a library, or perhaps it is C code.

Another good use of interfaces are to create generic subroutines and for overloading operators. For example, you could use an interface to overide the Fortran multiplication operator * to do the cross-product in the above example. So instead of using the CrossProduct function shown above you could just write c = a * b.

--
David

Re: Problem declaring function which outputs vector

By way of example, one of many modules I am using reads contains, inter alia the following:

      MODULE LOADING
C     --------------
C     jw 25-05-12 draft
C     jw 25-05-12 last amended

C     INTEGER  NLOAD  !Declared in module PARAMS (file: mPARAMS)
C     INTEGER  NLCAS  !Declared in module PARAMS (file: mPARAMS)
C     INTEGER  NCOMB  !Declared in module PARAMS (file: mPARAMS)
C     INTEGER  NLCASC !Declared in module PARAMS (file: mPARAMS)

      CHARACTER*48, ALLOCATABLE :: LCName(:), CombName(:)

      Type :: LOADS
         REAL             ::   Pos
         CHARACTER(36)    ::   Desc
         CHARACTER(1)     ::   Axes, LTyp, TDLshp
         REAL             ::   Fx,Fy,Mz
      END TYPE

      TYPE(LOADS),   ALLOCATABLE :: LOADCASE(:,:)

      Type :: COMBS
         INTEGER          ::   LCnum
         REAL             ::   Factr
      END TYPE

      TYPE(COMBS),   ALLOCATABLE :: LOADCOMB(:,:)

      END MODULE Loading

I have USE statements for this in all program files which use these arrays, including, of course, the program files in which the arrays are initially defined and then initialised.

Example:

      SUBROUTINE ALLOC2
C     -----------------
C     jw / 25-05-12  first drafted
C     jw / 05-06-12  last rev.
C
C     ALLOCATE dynamic array sizes for loads

      USE PARAMS
      USE LOADING

C     IF(NLOAD.LT.NNODE/2 .OR. NLOAD.GT.NNODE*NMEMB)
C    +   NLOAD=(NMEMB+NNODE)*2 !Default no. of load entries per loadcase

      IF(ALLOCATED(LOADCASE)) DEALLOCATE(LOADCASE)
      IF(ALLOCATED(LCName)) DEALLOCATE(LCName)

      ALLOCATE(LCName(NLCAS))
      ALLOCATE(LOADCASE(NLCAS,NLOAD))

      RETURN
      END

Is the first of these by any chance a form of an 'interface' module?
---
John

7 Reply by davidb (edited by davidb 2015-01-15 07:23:47)

Re: Problem declaring function which outputs vector

JohnWasilewski wrote:

Is the first of these by any chance a form of an 'interface' module?

Not really. It is a module which defines some derived types and arrays which are shared amongst subroutines which USE the module. There are no interface blocks in your module.

This is an example of a module that uses an interface to create a sort subroutine that can handle real and double precision arrays automatically (and a derived type). Note that private is used to disallow calls to sort_real, sort_double and sort_foo_t directly.

module sorting

   private :: sort_real, sort_double, sort_foo_t

   type foo_t
     real :: a
     real :: b
     real :: key
   end type foo_t

   interface sort
      module procedure sort_real
      module procedure sort_double
      module procedure sort_foo_t
   end interface sort
   
contains

   !** Insertion sort of array of double precision values
   subroutine sort_double(a)
      real(kind=8), intent(inout) :: a(:)
      integer :: i, j
      real(kind=8) :: tmp
      print *,'sort_double called'
      do i=1, size(a)
            tmp = a(i)
            do j=i-1, 1, -1
               if (tmp >= a(j)) exit
               a(j+1) = a(j)
            end do
            a(j+1) = tmp
         end do
   end subroutine sort_double
   
   !** Insertion sort of array of real values
   subroutine sort_real(a)
      real, intent(inout) :: a(:)
      integer :: i, j
      real :: tmp
      print *,'sort_real called'
      do i=1, size(a)
            tmp = a(i)
            do j=i-1, 1, -1
               if (tmp >= a(j)) exit
               a(j+1) = a(j)
            end do
            a(j+1) = tmp
         end do
   end subroutine sort_real
   
   !** Insertion sort of array of foo_type(s)
   subroutine sort_foo_t(a)
      type(foo_t), intent(inout) :: a(:)
      integer :: i, j
      type(foo_t) :: tmp
      print *,'sort_foo_t called'
      do i=1, size(a)
            tmp = a(i)
            do j=i-1, 1, -1
               if (tmp%key >= a(j)%key) exit
               a(j+1) = a(j)
            end do
            a(j+1) = tmp
         end do
   end subroutine sort_foo_t
   
end module

program main

   use sorting
   
   real :: x(6) = (/2.0, 4.0, 1.0, 0.0, 5.0, 3.0/)
   real(kind=8) :: y(6) = (/2.0d0, 4.0d0, 1.0d0, 0.0d0, 5.0d0, 3.0d0/)
   type(foo_t) :: z(6) = (/foo_t(0.0, 0.0, 2.0), &
                           foo_t(0.0, 0.0, 4.0), &
                           foo_t(0.0, 0.0, 1.0), &
                           foo_t(0.0, 0.0, 0.0), &
                           foo_t(0.0, 0.0, 5.0), &
                           foo_t(0.0, 0.0, 3.0)/)
                             
   ! Interface assures real version is called correctly.
   call sort(x)
   print *, x
   
   ! Interface assures double precision version is called correctly
   call sort(y)
   print *, y
   
   ! Interface assures foo_t version is called correctly
   call sort(z)
   print *, z%key

end program main
--
David