Cross Product

[St-Maxwell]

Hello, I tried to implement cross_product function as a component of linalg module. And a brief description is presented below. If we can have an agreement on the interface of cross_product, I will create a pull request.

Cross Product

Overview

Cross product is a binary operation on two vectors in space, denoted by . The result of is a vector that is perpendicular to both and . Let , , and being the orthonormal basis vectors in space, and can be represented by the following expressions


Then has a determinant form:

API

interface
    function cross_product(a, b) result(c)
        real, intent(in) :: a(:)
        real, intent(in) :: b(:)
        real :: c(3)
    end function
end interface

error stop when the size of a or b is not equal to 3.

Discussion

I think the vectors a and b should be passed as assumed-shape array. There is no need to explicitly define a and b as fixed-size array, i.e. real, intent(in) :: a(3), b(3). Compilers are not able to check whether the size is reasonable at compile-time since in many cases, allocatable arrays are used. Therefore, the assertion for array size is placed in the body of cross_product function, where error_stop is called. As a result, cross_product is no longer a pure function.

Based on the above discussion and consideration of simplicity, I chose an analogous way to Julia to implement cross_product.

Reference Implementations

• Numpy
• Julia
• C++ (Boost)

[ivan-pi]

There is a blog-post by @jacobwilliams on this function: http://degenerateconic.com/old-school/

I would be happy to drop the assertion and simply document the size of the arguments should be 3. Since it is a low-level function that might be called hundreds of times I think we don't want an error check.

[St-Maxwell]

@ivan-pi Thanks for sharing your idea. I agree that the assertion would lower the efficiency of cross_product. But I'm wondering if we could have a better solution, for example, a fypp preprocess statement that enables the assertion when debug and drops it for release.

[ivan-pi]

Having a preprocessor warning block would be one way, but I really think the cross-product should be a pure function to benefit from compiler inlining.

With the debug build in fpm or CMake we already have -fcheck=bounds for gfortran and -check bounds for ifort that should be able to capture if size(a) /= size(b) or if a and b are not allocated.

Do you have any suggestion for an interface to deal with the cross-product of a vector field, e.g. to solve the shallow water equations on a sphere:

integer, parameter :: n = 1000 ! number of grid points
real, allocatable :: x(:,:), u(:,:)

! Nodes on the sphere
allocate(x(3,n))

! Velocity vectors on the sphere surface
allocate(u(3,n))

Alternatively, I might decide to use a different memory layout:

! Array of Structures layout
allocate(u(n,3), x(n,3))

! or, AoS more explicitly
type :: vector_field_t
  integer :: n
  real, allocatable :: x(:), y(:), z(:)
end type

[ivan-pi]

I just wanted to add that requests for cross product have appeared previously on both comp.lang.fortran and StackExchange.