Computing the cross product of two vectors in Fortran 90

Question

I would like to compute the cross product of two vectors in Fortran 90. For example, in words, the cross product of (1, 2, 3) and (4, 5, 6) turns out to be (-3, 6, -3) in Cartesian coordinates. I wrote the following code (main program followed by function definition):

PROGRAM crosstest
  IMPLICIT NONE

  INTEGER, DIMENSION(3) :: m, n
  INTEGER, DIMENSION(3) :: cross
  INTEGER, DIMENSION(3) :: r

  m=(/1, 2, 3/)
  n=(/4, 5, 6/)
  r=cross(m,n)

END PROGRAM crosstest

FUNCTION cross(a, b)
  INTEGER, DIMENSION(3) :: cross
  INTEGER, DIMENSION(3), INTENT(IN) :: a, b

  cross(1) = a(2) * b(3) - a(3) * b(2)
  cross(2) = a(3) * b(1) - a(1) * b(3)
  cross(3) = a(1) * b(2) - a(2) * b(1)
END FUNCTION cross

But, I get an error message:

crosstest.f90:10.9:

  r=cross(m,n)
         1
Error: Rank mismatch in array reference at (1) (2/1)

where line 10 is r=cross(m,n). It seems that I must be specifying a dimension incorrectly. Here are a few ideas I have:

1. Perhaps the declaration of the function cross in the main program should be simply an integer variable, rather than a 1by3 integer array. So I tried deleting the , DIMENSION(3) in the INTEGER, DIMENSION(3) :: cross line in the main program. But I get an error message:

   crosstest.f90:10.4:
   
     r=cross(m,n)
       1
   Error: The reference to function 'cross' at (1) either needs an
   explicit INTERFACE or the rank is incorrect
   

   so this is even worse, probably.

2. Some (but not all) Fortran function examples on the web place an EXTERNAL statement after the function declaration in the main program. So I tried placing a line EXTERNAL cross after the declaration block in the main program. I get an error message:

   crosstest.f90:8.16:
   
     EXTERNAL cross
                   1
   Error: EXTERNAL attribute conflicts with DIMENSION attribute at (1)
   

   So this seems incorrect also.

3. Some (but not all) Fortran function examples on the web place a RETURN statement on the second-to-last line of the function definition. I tried this, but I get the original rank mismatch error:

   crosstest.f90:10.9:
   
     r=cross(m,n)
            1
   Error: Rank mismatch in array reference at (1) (2/1)
   

   So this does not fix the problem.

Can you please help me see my error?

Urgency

Answer 1

The best practice is to place your procedures (subroutines and functions) in a module and then "use" that module from your main program or other procedures. You don't need to "use" the module from other procedures of the same module. This will make the interface of the procedure explicit so that the calling program or procedure "knows" the characteristics of the arguments ... it allows the compiler to check for consistency between the arguments on both sides ... caller and callee .. this eliminates a lot of bugs.

Outside of the language standard, but in practice necessary: if you use one file, place the module before the main program that uses it. Otherwise the compiler will be unaware of it. so:

module my_subs

implicit none

contains

FUNCTION cross(a, b)
  INTEGER, DIMENSION(3) :: cross
  INTEGER, DIMENSION(3), INTENT(IN) :: a, b

  cross(1) = a(2) * b(3) - a(3) * b(2)
  cross(2) = a(3) * b(1) - a(1) * b(3)
  cross(3) = a(1) * b(2) - a(2) * b(1)
END FUNCTION cross

end module my_subs


PROGRAM crosstest
  use my_subs
  IMPLICIT NONE

  INTEGER, DIMENSION(3) :: m, n
  INTEGER, DIMENSION(3) :: r

  m= [ 1, 2, 3 ]
  n= [ 4, 5, 6 ]
  r=cross(m,n)
  write (*, *) r

END PROGRAM crosstest

Whitcomb

Answer 2

This is kind of a late answer, but since I stumbled upon this and there is no real explanation yet for why your error occurred, I figured I'd add an explanation for everybody else who stumbles upon this question:

In your program, you define an array called cross, which is of rank 1. Then you call the cross function you define further down. Since the cross function does not have an explicit interface (see M.S.B.'s answer), the compiler does not know about it at this point. What it does know about is the array you declared. If you write r = cross(m, n), the compiler thinks you want to access the element at position (m, n) of the array cross. Since this array is of rank 1, but you supplied two arguments, you get the error

rank mismatch in array reference at (1) (2/1)

which means that you supplied two coordinates when the compiler was expecting one.

Karp

Thanks. That makes sense. I edited the answer accordingly. – Karp

Answer 3

You can place the subroutines used in a program after a contains keyword within the program. This eliminates the need for creating a module or adding the interface definition.

PROGRAM crosstest
  IMPLICIT NONE

  INTEGER, DIMENSION(3) :: m, n
  INTEGER, DIMENSION(3) :: r

  m = (/1, 2, 3/)
  n = (/4, 5, 6/)
  r = cross(m,n)

  print *, r

CONTAINS

PURE FUNCTION cross(a, b)
  INTEGER, DIMENSION(3) :: cross
  INTEGER, DIMENSION(3), INTENT(IN) :: a, b

  cross(1) = a(2) * b(3) - a(3) * b(2)
  cross(2) = a(3) * b(1) - a(1) * b(3)
  cross(3) = a(1) * b(2) - a(2) * b(1)
END FUNCTION cross

END PROGRAM crosstest

Limousin

INTEGER, DIMENSION(3) :: cross This line should be deleted, otherwise compiling errors appear – Criticize

@YoujunHu - you are correct. Fixed. – Limousin

Answer 4

The complier interpretes the cross in your main program as an array, rather than a function name. Therefore cross in your main program has nothing to do with the function you defined. To let the compiler know that cross is a function name, you can just define in your main program:

INTEGER :: cross

Then, when the compiler see the line r=cross(m,n), the compiler know that this is not array indexing since cross is not defined as an array in the main program. Then a remaining possiblity is that this is a function call, and Fortran compiler ususally adopt this possiblity.

For your case, the compiler finds out that the cross function returns an array rather than a scalar. For this case, the compiler insists that you provide an explicit interface for cross. So this workaround does not work for your case. As other anwsers suggest, you need to provide an explicit interface by either

• containing the function definition in a module and use that module in your main program

or

• containing the function definition in your main program.

Fortran uses parentheses for both function call and array indexing. This design often causes confusion because a programmer can not figure out f(3) is a function call or array indexing without additional information.

Criticize

It is poor advice, the suggestion to declare cross as an integer in the main program. Which is a little odd, as the advice later in the answer, to wrap the function into a module or to contain it inside the main program, is good advice. Albeit it repeats the meat of the already accepted answer. – Unwatched

That is not advice, but reasoning why the compiler can not find the function and how to let the compiler find the function within the framework represented by the code in the question. And this is a common practice in Fortran77 when one wants to use an external function. So it is not odd. – Criticize