Determinant in Fortran95

Question

This code in fortran calculates the determinant of a nxn matrix using the laplacian formula (expansion by minors). I understand fully how this process works.

But could somebody give me an insight into how the following code operates over, say a given iteration, this section of the code contains the recursive function determinant(matrix) - assume some nxn matrix is read in and passed through and another function to call the cofactor. There are aspects of the code I understand but its the recursion that is confusing me profoundly. I've tried to run through step by step with a 3x3 matrix but to no avail.

! Expansion of determinants using Laplace formula
recursive function determinant(matrix) result(laplace_det)
real, dimension(:,:) :: matrix
integer :: msize(2), i, n
real :: laplace_det, det
real, dimension(:,:), allocatable :: cf

msize = shape(matrix) 
n = msize(1)          

if (n .eq. 1) then  
  det = matrix(1,1)
else
  det = 0    
  do i=1, n  
    allocate(cf(n-1, n-1))     
    cf = cofactor(matrix, i, 1)
    det = det + ((-1)**(i+1))* matrix(i,1) * determinant(cf)
    deallocate(cf)
  end do        
end if
laplace_det = det
end function determinant       

 function cofactor(matrix, mI, mJ)
  real, dimension(:,:) :: matrix
  integer :: mI, mJ
  integer :: msize(2), i, j, k, l, n
  real, dimension(:,:), allocatable :: cofactor
  msize = shape(matrix)
  n = msize(1)

  allocate(cofactor(n-1, n-1))
  l=0
  k = 1
  do i=1, n
   if (i .ne. mI) then
     l = 1
     do j=1, n
       if (j .ne. mJ) then
         cofactor(k,l) = matrix(i,j)
         l = l+ 1
       end if
     end do
     k = k+ 1
   end if
  end do
return
end function cofactor

The main section im struggling with is these two calls and the operation of the respective cofactor calculation.

cf = cofactor(matrix, i, 1)
det = det + ((-1)**(i+1))* matrix(i,1) * determinant(cf)

Any input for an explanation would be greatly appreciated (like i said an example of one iteration). This is my first post in stack-overflow as most of my question reside in mathstack (as you can probably tell by the mathematical nature of the question). Although I do have experience programming, the concept of recursion (especially in this example) is really boggling my mind.

If any edits are required please go ahead, im not familiar with the etiquette on stack overflow.

Answer 1

Let us suppose that we pass the following 3x3 matrix to determinant():

2 9 4
7 5 3
6 1 8

In the routine, the following two lines are executed iteratively for i = 1,2,3:

cf = cofactor(matrix, i, 1)
det = det + ((-1)**(i+1))* matrix(i,1) * determinant(cf)

which corresponds to the Laplace expansion with respect to the first column. More specifically, one passes the above 3x3 matrix to cofactor() to get a 2x2 sub-matrix by removing the i-th row and 1st column of the matrix. The obtained 2x2 sub-matrix (cf) is then passed to determinant() in the next line to calculate the co-factor corresponding to this sub-matrix. So, in this first iterations we are trying to calculate

Note here that the three determinants in the right-hand side are yet to be calculated by subsequent calls of determinant(). Let us consider one such subsequent call, e.g. for i=1. We are passing the following sub-matrix (stored in cf)

5 3
1 8

to determinant(). Then, the same procedure as described above is repeated again and independently of the Laplace expansion for the parent 3x3 matrix. That is, the determinant() now iterates over i=1,2 and tries to calculate

Note that the i in this subsequent call is different from the i of the previous call; they are all local variables living inside a particular call of a routine and are totally independent from each other. Also note that the index of dummy array argument (like matrix(:,:)) always start from 1 in Fortran (unless otherwise specified). This kind of operations are repeated until the size of the sub-matrix becomes 1.

But in practice, I believe that the easiest way to understand this kind of code is to write intermediate data and track the actual flow of data/routines. For example, we can insert a lot of print statements as

module mymod
    implicit none
contains

recursive function determinant(matrix) result(laplace_det)
    real :: matrix(:,:)
    integer :: i, n, p, q
    real :: laplace_det, det
    real, allocatable :: cf(:,:)

    n = size(matrix, 1) 

    !***** output *****   
    print "(a)", "Entering determinant() with matrix = "
    do p = 1, n
        print "(4x,100(f3.1,x))", ( matrix( p, q ), q=1,n )
    enddo

    if (n == 1) then  
        det = matrix(1,1)
    else
        det = 0
        do i = 1, n  
            allocate( cf(n-1, n-1) )
            cf = cofactor( matrix, i, 1 )

            !***** output *****
            print "(4x,a,i0,a,i0,a)", "Getting a ", &
                    n-1, "-by-", n-1, " sub-matrix from cofactor():"
            do p = 1, n-1
                print "(8x, 100(f3.1,x))", ( cf( p, q ), q=1,n-1 )
            enddo

            print "(4x,a)", "and passing it to determinant()."

            det = det + ((-1)**(i+1))* matrix( i, 1 ) * determinant( cf )
            deallocate(cf)
        end do
    end if

    laplace_det = det

    !***** output *****
    print *, "  ---> Returning det = ", det
end function

function cofactor(matrix, mI, mJ)
    .... (same as the original code)
end function

end module

program main
    use mymod
    implicit none
    real :: a(3,3), det

    a( 1, : ) = [ 2.0, 9.0, 4.0 ]
    a( 2, : ) = [ 7.0, 5.0, 3.0 ]
    a( 3, : ) = [ 6.0, 1.0, 8.0 ]

    det = determinant( a )
    print "(a, es30.20)", "Final det = ", det
end program

then the output clearly shows how the data are processed:

Entering determinant() with matrix = 
    2.0 9.0 4.0
    7.0 5.0 3.0
    6.0 1.0 8.0
    Getting a 2-by-2 sub-matrix from cofactor():
        5.0 3.0
        1.0 8.0
    and passing it to determinant().
Entering determinant() with matrix = 
    5.0 3.0
    1.0 8.0
    Getting a 1-by-1 sub-matrix from cofactor():
        8.0
    and passing it to determinant().
Entering determinant() with matrix = 
    8.0
   ---> Returning det =    8.0000000    
    Getting a 1-by-1 sub-matrix from cofactor():
        3.0
    and passing it to determinant().
Entering determinant() with matrix = 
    3.0
   ---> Returning det =    3.0000000    
   ---> Returning det =    37.000000    
    Getting a 2-by-2 sub-matrix from cofactor():
        9.0 4.0
        1.0 8.0
    and passing it to determinant().
Entering determinant() with matrix = 
    9.0 4.0
    1.0 8.0
    Getting a 1-by-1 sub-matrix from cofactor():
        8.0
    and passing it to determinant().
Entering determinant() with matrix = 
    8.0
   ---> Returning det =    8.0000000    
    Getting a 1-by-1 sub-matrix from cofactor():
        4.0
    and passing it to determinant().
Entering determinant() with matrix = 
    4.0
   ---> Returning det =    4.0000000    
   ---> Returning det =    68.000000    
    Getting a 2-by-2 sub-matrix from cofactor():
        9.0 4.0
        5.0 3.0
    and passing it to determinant().
Entering determinant() with matrix = 
    9.0 4.0
    5.0 3.0
    Getting a 1-by-1 sub-matrix from cofactor():
        3.0
    and passing it to determinant().
Entering determinant() with matrix = 
    3.0
   ---> Returning det =    3.0000000    
    Getting a 1-by-1 sub-matrix from cofactor():
        4.0
    and passing it to determinant().
Entering determinant() with matrix = 
    4.0
   ---> Returning det =    4.0000000    
   ---> Returning det =    7.0000000    
   ---> Returning det =   -360.00000    
Final det =    -3.60000000000000000000E+02

You can insert more print lines until the whole mechanism becomes clear.

---

BTW, the code in the Rossetta page seems much simpler than the OP's code by creating a sub-matrix directly as a local array. The simplified version of the code reads

recursive function det_rosetta( mat, n ) result( accum )
    integer :: n
    real    :: mat(n, n)
    real    :: submat(n-1, n-1), accum
    integer :: i, sgn

    if ( n == 1 ) then
        accum = mat(1,1)
    else
        accum = 0.0
        sgn = 1
        do i = 1, n
            submat( 1:n-1, 1:i-1 ) = mat( 2:n, 1:i-1 )
            submat( 1:n-1, i:n-1 ) = mat( 2:n, i+1:n )

            accum = accum + sgn * mat(1, i) * det_rosetta( submat, n-1 )
            sgn = - sgn
        enddo
    endif
end function

Note that the Laplace expansion is made along the first row, and that the submat is assigned using array sections. The assignment can also be written simply as

submat( :, :i-1 ) = mat( 2:, :i-1 )
submat( :, i: )   = mat( 2:, i+1: )

where the upper and lower bounds of the array sections are omitted (then, the declared values of upper and lower bounds are used by default). The latter form is used in the Rosetta page.