Cross Product In Matrix Form

Motivation

A proposal to include a cross_product_matrix() function in stdlib_linalg module that takes a 3×1 vector and returns the 3×3 skew-symmetric cross product operator matrix which is used in mechanics a lot to convert the cross product into a linear algebra operator.

An illustrative example of this function would be:

    pure function cross_product_matrix(a) result(c)
    real(dp), intent(in) ::a(3)
    reaL(dp) :: c(3,3)
    real(dp) :: x,y,z

    x = a(1)
    y = a(2)
    z = a(3)
    ! define column-by-column (negative of row-major display)
    c = reshape( &
        [ 0._dp,  z,                -y, &
                  -z, 0._dp,         x, &
                   y,        -x, 0._dp], [3,3])
    end function cross_product_matrix

This allows a vector cross producrt c=a×b to be described using linear algebra as the matrix-vector multiplication c=A*b where A=a× is the cross product matrix.

The use of this function would be in mechanics. For example building a 3×3 rotation matrix from an axis vector and an angle scalar

    pure function rotation_matrix(axis, angle) result(R)
    real(dp) :: R(3,3)
    real(dp), intent(in) :: angle, axis(3)
    real(dp) :: ax(3,3)

        ! Form an identiy matrix
        R = 0._dp
        R(1,1) = 1._dp
        R(2,2) = 1._do
        R(3,3) = 1._dp

        ! Form the cross product matrix
        ax = vector_cross_matrix(axis)

        ! Rodrigues Rotation Formula
        R = R + sin(angle)*ax + (1._dp - cos(angle))*matmul(ax, ax)

    end function rotation_matrix

or finding the mass moment of inertia tensor of a body based on the parallel axis theorem

    pure function parallel_axis_theorem(I, m, r) result(Ip)
    real(dp) :: Ip(3,3)
    real(dp), intent(in) :: I(3,3), m, r(3)
    real(dp) :: rx(3,3)

        ! Form the cross product matrix
        rx = vector_cross_matrix(r)

        ! Parallel Axis Theorem
        Ip = I - m*(matmul(rx, rx))
    end function parallel_axis_theorem

Prior Art

Existing cross product function cross_product() takes two vectors and returns their cross product vector as a result.

stdlib/src/stdlib_linalg_cross_product.fypp

Line 4 in 2bdc50e

and the linear algebra module definition file:

stdlib/src/stdlib_linalg.fypp

Line 4 in 2bdc50e

#:set RHS_SUFFIX = ["one","many"]

Additional Information

No response