Déterminant dans Fortran95

Question

Ce code dans fortran calcule le déterminant d'une matrice nxn en utilisant la formule laplacienne (expansion par des mineurs). Je comprends parfaitement comment ce processus fonctionne. Mais quelqu'un pourrait-il me donner un aperçu du fonctionnement du code suivant, disons une itération donnée, cette section du code contient le déterminant de la fonction récursive (matrice) - supposons qu'une matrice nxn est lue et transmise et une autre fonction pour appeler le cofacteur. Il y a des aspects du code que je comprends mais c'est la récursivité qui me trouble profondément. J'ai essayé de passer à travers étape par étape avec une matrice 3x3 mais en vain.

! 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 

La section principale avec laquelle nous luttons est ces deux appels et le fonctionnement du calcul du cofacteur respectif. Toute entrée pour une explication serait grandement appréciée (comme j'ai dit un exemple d'une seule itération). Ceci est mon premier article en pile-débordement comme la plupart de ma question réside dans mathstack (comme vous pouvez probablement le dire par la nature mathématique de la question). Même si j'ai de l'expérience en programmation, le concept de récursivité (surtout dans cet exemple) est vraiment ahurissant.

Si des modifications sont nécessaires, veuillez aller de l'avant, je ne suis pas familier avec l'étiquette sur le débordement de la pile.

Answer 1

Supposons que nous passons la matrice 3x3 suivante determinant():

2 9 4 
7 5 3 
6 1 8 

Dans la routine, les deux lignes suivantes sont exécutées itérativement pour i = 1,2,3:

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

qui correspond à la Laplace expansion par rapport à la première colonne. Plus précisément, on passe les 3x3 matrix ci-dessus à cofactor() pour obtenir une sous-matrice 2x2 en supprimant la 0e ligne i et la 1ère colonne du matrix. La sous-matrice 2x2 obtenue (cf) est ensuite passée à determinant() dans la ligne suivante pour calculer le cofacteur correspondant à cette sous-matrice. Ainsi, dans ce premier itérations que nous essayons de calculer

Notez ici que les trois déterminants du côté droit ne sont pas encore calculées par des appels ultérieurs de déterminant(). Considérons un tel appel ultérieur, par ex. pour i=1. Nous passons la sous-matrice suivante (stockée dans cf)

5 3 
1 8 

à determinant(). Ensuite, la même procédure que celle décrite ci-dessus est répétée à nouveau et indépendamment de l'expansion de Laplace pour la matrice mère 3x3. Autrement dit, le déterminant() maintenant itère sur i=1,2 et essaie de calculer

Notez que le i dans cet appel ultérieur est différent du i de l'appel précédent; ce sont toutes des variables locales qui vivent à l'intérieur d'un appel particulier d'une routine et sont totalement indépendantes les unes des autres.Notez également que l'index de l'argument du tableau factice (comme matrix(:,:)) commence toujours par 1 dans Fortran (sauf indication contraire). Ce type d'opérations est répété jusqu'à ce que la taille de la sous-matrice devienne 1. Mais en pratique, je crois que la façon la plus simple de comprendre ce type de code est d'écrire des données intermédiaires et de suivre le flux réel des données/routines. Par exemple, nous pouvons insérer beaucoup de print déclarations comme

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 

alors la sortie montre clairement comment les données sont traitées:

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 

Vous pouvez insérer plusieurs lignes d'impression jusqu'à ce que tout le mécanisme est clair.

---

BTW, le code the Rossetta page semble beaucoup plus simple que le code de l'OP en créant une sous-matrice directement en tant que tableau local. La version simplifiée du code lit

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 

Notez que l'expansion de Laplace est faite le long de la première ligne, et que la submat est affectée à l'aide des sections de tableau. L'affectation peut aussi être écrit simplement comme

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

où les limites supérieures et inférieures des sections de tableau sont omis (alors, les valeurs déclarées de limites supérieures et inférieures sont utilisées par défaut). Ce dernier formulaire est utilisé dans la page Rosetta.