让我们假设我们将以下 3x3 矩阵传递给determinant()：
2 9 4
7 5 3
6 1 8
在例程中，迭代执行以下两行代码i = 1,2,3：
cf = cofactor(matrix, i, 1)
det = det + ((-1)**(i+1))* matrix(i,1) * determinant(cf)
这对应于第一列的拉普拉斯展开。更具体地说，通过删除上述 3x3的第 - 行和第 1 列matrix来cofactor()获得 2x2 子矩阵。然后将获得的 2x2 子矩阵 ( ) 传递到下一行，以计算对应于该子矩阵的 co-factor。所以，在第一次迭代中，我们试图计算imatrixcfdeterminant()
请注意，右侧的三个行列式尚未通过随后的行列式 () 调用来计算。让我们考虑一个这样的后续调用，例如 for i=1。我们正在传递以下子矩阵（存储在cf）
5 3
1 8
到determinant(). 然后，与父 3x3 矩阵的拉普拉斯展开无关，再次重复上述相同过程。也就是说，determinant() 现在迭代i=1,2并尝试计算
请注意，i此后续调用中i的 与前一次调用中的不同；它们都是存在于例程的特定调用中的局部变量，并且彼此完全独立。另请注意，虚拟数组参数（如matrix(:,:)）的索引始终从1Fortran 开始（除非另有说明）。重复这种操作，直到子矩阵的大小变为1。
但在实践中，我认为理解这类代码的最简单方法是编写中间数据并跟踪数据/例程的实际流动。例如，我们可以插入很多print语句为
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
然后输出清楚地显示了数据是如何处理的：
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
您可以插入更多打印行，直到整个机制变得清晰。
顺便说一句，通过直接将子矩阵创建为本地数组，Rossetta 页面中的代码似乎比 OP 的代码简单得多。代码的简化版本如下
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
请注意，拉普拉斯展开是沿第一行进行的，并且submat是使用数组部分分配的。作业也可以简单地写成
submat( :, :i-1 ) = mat( 2:, :i-1 )
submat( :, i: ) = mat( 2:, i+1: )
其中数组部分的上限和下限被省略（然后，默认使用上限和下限的声明值）。后一种形式用于 Rosetta 页面。