Problem Statement

Using the most straightfoward algorithm (which we assume here), computing the product of two matrices of dimensions (n1,n2) and (n2,n3) requires n1n2n3 FMA operations. The number of operations required to compute the product of matrices A1, A2... An depends on the order of matrix multiplications, hence on where parens are put. Remember that the matrix product is associative, but not commutative, hence only the parens can be moved. For instance, with four matrices, one can compute A(B(CD)), A((BC)D), (AB)(CD), (A(BC))D, (AB)C)D. The number of different ways to put the parens is a Catalan number, and grows exponentially with the number of factors. Here is an example of computation of the total cost, for matrices A(5,6), B(6,3), C(3,1): In this case, computing (AB)C requires more than twice as many operations as A(BC). The difference can be much more dramatic in real cases. Write a function which, given a list of the successive dimensions of matrices A1, A2... An, of arbitrary length, returns the optimal way to compute the matrix product, and the total cost. Any sensible way to describe the optimal solution is accepted. The input list does not duplicate shared dimensions: for the previous example of matrices A,B,C, one will only pass the list [5,6,3,1] (and not [5,6,6,3,3,1]) to mean the matrix dimensions are respectively (5,6), (6,3) and (3,1). Hence, a product of n matrices is represented by a list of n+1 dimensions. Try this function on the following two lists: To solve the task, it's possible, but not required, to write a function that enumerates all possible ways to parenthesize the product. This is not optimal because of the many duplicated computations, and this task is a classic application of dynamic programming. See also Matrix chain multiplication on Wikipedia.

Let's start with the solution:

module optim_mod
    implicit none
contains
    subroutine optim(a)
        implicit none
        integer :: a(:), n, i, j, k
        integer, allocatable :: u(:, :)
        integer(8) :: c
        integer(8), allocatable :: v(:, :)

        n = ubound(a, 1) - 1
        allocate (u(n, n), v(n, n))
        v = huge(v)
        u(:, 1) = -1
        v(:, 1) = 0
        do j = 2, n
            do i = 1, n - j + 1
                do k = 1, j - 1
                    c = v(i, k) + v(i + k, j - k) + int(a(i), 8) * int(a(i + k), 8) * int(a(i + j), 8)
                    if (c < v(i, j)) then
                        u(i, j) = k
                        v(i, j) = c
                    end if
                end do
            end do
        end do
        write (*, "(I0,' ')", advance="no") v(1, n)
        call aux(1, n)
        print *
        deallocate (u, v)
    contains
        recursive subroutine aux(i, j)
            integer :: i, j, k

            k = u(i, j)
            if (k < 0) then
                write (*, "(I0)", advance="no") i
            else
                write (*, "('(')", advance="no")
                call aux(i, k)
                write (*, "('*')", advance="no")
                call aux(i + k, j - k)
                write (*, "(')')", advance="no")
            end if
        end subroutine
    end subroutine
end module

program matmulchain
    use optim_mod
    implicit none

    call optim([5, 6, 3, 1])
    call optim([1, 5, 25, 30, 100, 70, 2, 1, 100, 250, 1, 1000, 2])
    call optim([1000, 1, 500, 12, 1, 700, 2500, 3, 2, 5, 14, 10])
end program