Определитель в Fortran95

Предположим, что мы передаем следующую матрицу 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 matrix в cofactor(), чтобы получить подматрицу 2x2, удалив i-ю строку и 1-й столбец matrix. Полученная подматрица 2x2 (cf) затем передается в determinant() в следующей строке для вычисления кофактора, соответствующего этой подматрице. Итак, в этих первых итерациях мы пытаемся вычислить

Обратите внимание, что три определителя в правой части еще предстоит вычислить последующими вызовами determinant(). Рассмотрим один такой последующий вызов, например. для i=1. Мы передаем следующую подматрицу (хранится в cf)

5 3
1 8

до determinant(). Затем та же процедура, что описана выше, повторяется снова и независимо от разложения Лапласа для исходной матрицы 3x3. То есть определитель() теперь перебирает i=1,2 и пытается вычислить

Обратите внимание, что i в этом последующем вызове отличается от i предыдущего вызова; все они являются локальными переменными, живущими внутри конкретного вызова подпрограммы и полностью независимыми друг от друга. Также обратите внимание, что индекс аргумента фиктивного массива (например, matrix(:,:)) всегда начинается с 1 в Фортране (если не указано иное). Такого рода операции повторяются до тех пор, пока размер подматрицы не станет 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.