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.
