Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

Register Log in
  • Congratulations strongm on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Determinant using Laplace expansion

Status
Not open for further replies.
freddyn

freddyn

Technical User
Jan 22, 2013
3
DE
Hello!
Do somebody know where the code for calculating a determinant by Laplace theorem (in Fortran) can be found? Or maybe somebody could help with this issue.
Thank you!
 
To answer your question I tired this and it seems to work:

determinant.f95
Code: 
[COLOR=#a020f0]module[/color] matrix_functions
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#a020f0]contains[/color]
  [COLOR=#a020f0]subroutine[/color] print_matrix(matrix)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: matrix
    [COLOR=#2e8b57][b]integer[/b][/color] :: msize([COLOR=#ff00ff]2[/color]), i, j, m, n
    msize [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]shape[/color](matrix)
    m [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]1[/color])
    n [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]2[/color])
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], m
      [COLOR=#804040][b]do[/b][/color] j[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]'(f8.2)'[/color],[COLOR=#804040][b]advance[/b][/color] [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]'no'[/color]) matrix(i,j)
      [COLOR=#804040][b]end do[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])
    [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])
  [COLOR=#a020f0]end subroutine[/color] print_matrix 

  [COLOR=#a020f0]function[/color] cofactor(matrix, mI, mJ)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: matrix
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: mI, mJ
    [COLOR=#2e8b57][b]integer[/b][/color] :: msize([COLOR=#ff00ff]2[/color]), i, j, k, l, n
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]) :: cofactor
    msize [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]shape[/color](matrix)
    [COLOR=#804040][b]if[/b][/color] (msize([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b].ne.[/b][/color] msize([COLOR=#ff00ff]2[/color])) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Error in Cofactor: Input must be a square Matrix!'[/color]
      [COLOR=#804040][b]return[/b][/color]
    [COLOR=#804040][b]end if[/b][/color]
    n [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]1[/color])

    k [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
      [COLOR=#804040][b]if[/b][/color] (i [COLOR=#804040][b].ne.[/b][/color] mI) [COLOR=#804040][b]then[/b][/color]
         l [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
         [COLOR=#804040][b]do[/b][/color] j[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
           [COLOR=#804040][b]if[/b][/color] (j [COLOR=#804040][b].ne.[/b][/color] mJ) [COLOR=#804040][b]then[/b][/color]
             cofactor(k,l) [COLOR=#804040][b]=[/b][/color] matrix(i,j)
             l [COLOR=#804040][b]=[/b][/color] l[COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
           [COLOR=#804040][b]end if[/b][/color]
         [COLOR=#804040][b]end do[/b][/color]
         k [COLOR=#804040][b]=[/b][/color] k[COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
      [COLOR=#804040][b]end if[/b][/color]
    [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]return[/b][/color]
  [COLOR=#a020f0]end function[/color] cofactor

  [COLOR=#0000ff]! Expansion of determinants using Laplace formula[/color]
  [COLOR=#a020f0]recursive[/color] [COLOR=#a020f0]function[/color] determinant(matrix) [COLOR=#a020f0]result[/color](laplace_det)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: matrix
    [COLOR=#2e8b57][b]integer[/b][/color] :: msize([COLOR=#ff00ff]2[/color]), i, n, sgn
[COLOR=#2e8b57][b]    real[/b][/color] :: laplace_det, det
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]allocatable[/b][/color] :: cf

    msize [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]shape[/color](matrix)
    [COLOR=#804040][b]if[/b][/color] (msize([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b].ne.[/b][/color] msize([COLOR=#ff00ff]2[/color])) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Error in Determinant: Input must be a square Matrix!'[/color]
      [COLOR=#804040][b]return[/b][/color]
    [COLOR=#804040][b]end if[/b][/color]
    n [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]1[/color])
    
    [COLOR=#804040][b]if[/b][/color] (n [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#0000ff]! determinant of matrix 1 x 1, i.e. number[/color]
      det [COLOR=#804040][b]=[/b][/color] matrix([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]1[/color])
    [COLOR=#804040][b]else[/b][/color]
      det [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
      [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
        sgn [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]merge[/color]([COLOR=#ff00ff]1[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#008080]mod[/color](i,[COLOR=#ff00ff]2[/color]) [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]1[/color])
        [COLOR=#804040][b]allocate[/b][/color](cf(n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))
        cf [COLOR=#804040][b]=[/b][/color] cofactor(matrix, i, [COLOR=#ff00ff]1[/color])
        det [COLOR=#804040][b]=[/b][/color] det [COLOR=#804040][b]+[/b][/color] sgn [COLOR=#804040][b]*[/b][/color] matrix(i,[COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]*[/b][/color] determinant(cf)
        [COLOR=#804040][b]deallocate[/b][/color](cf)
      [COLOR=#804040][b]end do[/b][/color]        
    [COLOR=#804040][b]end if[/b][/color]
    [COLOR=#0000ff]![/color]
    laplace_det [COLOR=#804040][b]=[/b][/color] det
  [COLOR=#a020f0]end function[/color] determinant
[COLOR=#a020f0]end module[/color] matrix_functions


[COLOR=#a020f0]program[/color] main
  [COLOR=#a020f0]use[/color] matrix_functions

[COLOR=#2e8b57][b]  real[/b][/color] :: A1([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]1[/color]), A2([COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]2[/color]), A3([COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]3[/color]), A4([COLOR=#ff00ff]4[/color],[COLOR=#ff00ff]4[/color])
[COLOR=#2e8b57][b]  real[/b][/color] :: det 

  [COLOR=#0000ff]! matrix A1[/color]
  A1([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]1[/color])[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]5[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A1:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#a020f0]call[/color] print_matrix(A1)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A1:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A1)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A2[/color]
  A2([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]7[/color], [COLOR=#ff00ff]6[/color][COLOR=#804040][b]/[/b][/color])
  A2([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]9[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A2:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#a020f0]call[/color] print_matrix(A2)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A2:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A2)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A3[/color]
  A3([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]2[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A3([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]1[/color],  [COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A3([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]0[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A3:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#a020f0]call[/color] print_matrix(A3)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A3:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A3)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A4[/color]
  A4([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]2[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  A4([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]0[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color][COLOR=#804040][b]/[/b][/color])
  A4([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]6[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A4([COLOR=#ff00ff]4[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]5[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]2[/color],  [COLOR=#ff00ff]0[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A4:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#a020f0]call[/color] print_matrix(A4)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A4:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A4)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])  
[COLOR=#a020f0]end program[/color] main

Output:
Code: 
$ gfortran determinant.f95 -o determinant

$ determinant
 Matrix A1:
 ----------
    5.00

 Determinant of A1:
 -----------------
   5.0000000    

 Matrix A2:
 ----------
    7.00    6.00
    3.00    9.00

 Determinant of A2:
 -----------------
   45.000000    

 Matrix A3:
 ----------
   -2.00    2.00   -3.00
   -1.00    1.00    3.00
    2.00    0.00   -1.00

 Determinant of A3:
 -----------------
   18.000000    

 Matrix A4:
 ----------
    3.00    0.00    2.00   -1.00
    1.00    2.00    0.00   -2.00
    4.00    0.00    6.00   -3.00
    5.00    0.00    2.00    0.00

 Determinant of A4:
 -----------------
   20.000000
 
As mentioned here
this construct
Code: 
  function cofactor(matrix, mI, mJ)
    ...
    real, dimension(n-1, n-1) :: cofactor
    msize = shape(matrix)
    ...
    n = msize(1)
doesn't compile with newer gfortran version.
I tried it with g95 and got the similar error too:
Code: 
In file determinant.f95:23

    real, dimension(n-1, n-1) :: cofactor
                    1
Error: Variable 'n' cannot appear in restricted expression at (1)
It seems that it compiled only with my older gfortran version 4.4.0
 
The corrected code:
Code: 
[COLOR=#a020f0]module[/color] matrix_functions
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#a020f0]contains[/color]
  [COLOR=#a020f0]subroutine[/color] print_matrix(matrix)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: matrix
    [COLOR=#2e8b57][b]integer[/b][/color] :: msize([COLOR=#ff00ff]2[/color]), i, j, m, n
    msize [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]shape[/color](matrix)
    m [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]1[/color])
    n [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]2[/color])
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], m
      [COLOR=#804040][b]do[/b][/color] j[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
        [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#ff00ff]'(f8.2)'[/color],[COLOR=#804040][b]advance[/b][/color] [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]'no'[/color]) matrix(i,j)
      [COLOR=#804040][b]end do[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])
    [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])
  [COLOR=#a020f0]end subroutine[/color] print_matrix 

  [COLOR=#a020f0]function[/color] cofactor(matrix, mI, mJ)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: matrix
    [COLOR=#2e8b57][b]integer[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: mI, mJ
    [COLOR=#2e8b57][b]integer[/b][/color] :: msize([COLOR=#ff00ff]2[/color]), i, j, k, l, n
    [COLOR=#0000ff]! this doesn't work with with g95 and newer gfortran version: [/color]
    [COLOR=#0000ff]!real, dimension(n-1, n-1) :: cofactor[/color]
    [COLOR=#0000ff]! but this works with g95 and newer gfortran version:[/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]allocatable[/b][/color] :: cofactor
    msize [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]shape[/color](matrix)
    [COLOR=#804040][b]if[/b][/color] (msize([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b].ne.[/b][/color] msize([COLOR=#ff00ff]2[/color])) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Error in Cofactor: Input must be a square Matrix!'[/color]
      [COLOR=#804040][b]return[/b][/color]
    [COLOR=#804040][b]end if[/b][/color]
    n [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]1[/color])
    [COLOR=#0000ff]![/color]
    [COLOR=#804040][b]allocate[/b][/color](cofactor(n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))

    k [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
      [COLOR=#804040][b]if[/b][/color] (i [COLOR=#804040][b].ne.[/b][/color] mI) [COLOR=#804040][b]then[/b][/color]
         l [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
         [COLOR=#804040][b]do[/b][/color] j[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
           [COLOR=#804040][b]if[/b][/color] (j [COLOR=#804040][b].ne.[/b][/color] mJ) [COLOR=#804040][b]then[/b][/color]
             cofactor(k,l) [COLOR=#804040][b]=[/b][/color] matrix(i,j)
             l [COLOR=#804040][b]=[/b][/color] l[COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
           [COLOR=#804040][b]end if[/b][/color]
         [COLOR=#804040][b]end do[/b][/color]
         k [COLOR=#804040][b]=[/b][/color] k[COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
      [COLOR=#804040][b]end if[/b][/color]
    [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]return[/b][/color]
  [COLOR=#a020f0]end function[/color] cofactor

  [COLOR=#0000ff]! Expansion of determinants using Laplace formula[/color]
  [COLOR=#a020f0]recursive[/color] [COLOR=#a020f0]function[/color] determinant(matrix) [COLOR=#a020f0]result[/color](laplace_det)
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: matrix
    [COLOR=#2e8b57][b]integer[/b][/color] :: msize([COLOR=#ff00ff]2[/color]), i, n, sgn
[COLOR=#2e8b57][b]    real[/b][/color] :: laplace_det, det
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color](:,:), [COLOR=#2e8b57][b]allocatable[/b][/color] :: cf

    msize [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]shape[/color](matrix)
    [COLOR=#804040][b]if[/b][/color] (msize([COLOR=#ff00ff]1[/color]) [COLOR=#804040][b].ne.[/b][/color] msize([COLOR=#ff00ff]2[/color])) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Error in Determinant: Input must be a square Matrix!'[/color]
      [COLOR=#804040][b]return[/b][/color]
    [COLOR=#804040][b]end if[/b][/color]
    n [COLOR=#804040][b]=[/b][/color] msize([COLOR=#ff00ff]1[/color])
    
    [COLOR=#804040][b]if[/b][/color] (n [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#0000ff]! determinant of matrix 1 x 1, i.e. number[/color]
      det [COLOR=#804040][b]=[/b][/color] matrix([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]1[/color])
    [COLOR=#804040][b]else[/b][/color]
      det [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
      [COLOR=#804040][b]do[/b][/color] i[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], n
        sgn [COLOR=#804040][b]=[/b][/color] [COLOR=#008080]merge[/color]([COLOR=#ff00ff]1[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#008080]mod[/color](i,[COLOR=#ff00ff]2[/color]) [COLOR=#804040][b].eq.[/b][/color] [COLOR=#ff00ff]1[/color])
        [COLOR=#804040][b]allocate[/b][/color](cf(n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color], n[COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color]))
        cf [COLOR=#804040][b]=[/b][/color] cofactor(matrix, i, [COLOR=#ff00ff]1[/color])
        det [COLOR=#804040][b]=[/b][/color] det [COLOR=#804040][b]+[/b][/color] sgn [COLOR=#804040][b]*[/b][/color] matrix(i,[COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]*[/b][/color] determinant(cf)
        [COLOR=#804040][b]deallocate[/b][/color](cf)
      [COLOR=#804040][b]end do[/b][/color]        
    [COLOR=#804040][b]end if[/b][/color]
    [COLOR=#0000ff]![/color]
    laplace_det [COLOR=#804040][b]=[/b][/color] det
  [COLOR=#a020f0]end function[/color] determinant
[COLOR=#a020f0]end module[/color] matrix_functions


[COLOR=#a020f0]program[/color] main
  [COLOR=#a020f0]use[/color] matrix_functions

[COLOR=#2e8b57][b]  real[/b][/color] :: A1([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]1[/color]), A2([COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]2[/color]), A3([COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]3[/color]), A4([COLOR=#ff00ff]4[/color],[COLOR=#ff00ff]4[/color]), A5([COLOR=#ff00ff]5[/color],[COLOR=#ff00ff]5[/color])
[COLOR=#2e8b57][b]  real[/b][/color] :: det 

  [COLOR=#0000ff]! matrix A1[/color]
  A1([COLOR=#ff00ff]1[/color],[COLOR=#ff00ff]1[/color])[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]5[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A1:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#008080]call[/color] print_matrix(A1)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A1:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A1)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A2[/color]
  A2([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]7[/color], [COLOR=#ff00ff]6[/color][COLOR=#804040][b]/[/b][/color])
  A2([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]9[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A2:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#008080]call[/color] print_matrix(A2)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A2:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A2)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A3[/color]
  A3([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]2[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A3([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/-[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]1[/color],  [COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A3([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]0[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A3:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#008080]call[/color] print_matrix(A3)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A3:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A3)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A4[/color]
  A4([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]2[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  A4([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]0[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]2[/color][COLOR=#804040][b]/[/b][/color])
  A4([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]6[/color], [COLOR=#804040][b]-[/b][/color][COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A4([COLOR=#ff00ff]4[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]5[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]2[/color],  [COLOR=#ff00ff]0[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A4:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#008080]call[/color] print_matrix(A4)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A4:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A4)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])

  [COLOR=#0000ff]! matrix A5[/color]
  A5([COLOR=#ff00ff]1[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]5[/color], [COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]2[/color][COLOR=#804040][b]/[/b][/color])
  A5([COLOR=#ff00ff]2[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]1[/color], [COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]2[/color][COLOR=#804040][b]/[/b][/color])
  A5([COLOR=#ff00ff]3[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]6[/color], [COLOR=#ff00ff]2[/color], [COLOR=#ff00ff]3[/color][COLOR=#804040][b]/[/b][/color])
  A5([COLOR=#ff00ff]4[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]4[/color], [COLOR=#ff00ff]3[/color], [COLOR=#ff00ff]1[/color][COLOR=#804040][b]/[/b][/color])
  A5([COLOR=#ff00ff]5[/color],:)[COLOR=#804040][b]=[/b][/color]([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]0[/color], [COLOR=#ff00ff]2[/color][COLOR=#804040][b]/[/b][/color])
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Matrix A5:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'----------'[/color]
  [COLOR=#008080]call[/color] print_matrix(A5)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Determinant of A5:'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'-----------------'[/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) determinant(A5)
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color])    
[COLOR=#a020f0]end program[/color] main

Now it compiles fine with newer version of gfortran and g95:
Code: 
$ g95 determinant.f95 -o determinant

$ determinant
 Matrix A1:
 ----------
    5.00

 Determinant of A1:
 -----------------
 5.

 Matrix A2:
 ----------
    7.00    6.00
    3.00    9.00

 Determinant of A2:
 -----------------
 45.

 Matrix A3:
 ----------
   -2.00    2.00   -3.00
   -1.00    1.00    3.00
    2.00    0.00   -1.00

 Determinant of A3:
 -----------------
 18.

 Matrix A4:
 ----------
    3.00    0.00    2.00   -1.00
    1.00    2.00    0.00   -2.00
    4.00    0.00    6.00   -3.00
    5.00    0.00    2.00    0.00

 Determinant of A4:
 -----------------
 20.

 Matrix A5:
 ----------
    5.00    2.00    0.00    0.00    2.00
    0.00    1.00    4.00    3.00    2.00
    0.00    0.00    6.00    2.00    3.00
    0.00    0.00    4.00    3.00    1.00
    0.00    0.00    0.00    0.00    2.00

 Determinant of A5:
 -----------------
 100.
 
Thank you! Great job!
 
But the algorithm isn't very optimal, because the number of recursive calls seems to be n!
To reduce the number of recursions a little bit, you can compute the cases for N=2 and N=3 usning the known formulas without recursion - i.e. for N=2: m11*m22 - m12*m21 and for N=3: Rule of Sarus
 
yes thats a nice idea
 
Status
Not open for further replies.

Similar threads

R
  • Locked
  • Question
How do I print a real number using GNU gfortran
Replies
3
Views
423
rdflg
R
AlanJT
  • Locked
  • Question
I have developed a program using th
Replies
4
Views
157
AlanJT
AlanJT
alex422
  • Locked
  • Question
Floating-point exception using recursive function
Replies
10
Views
223
alex422
alex422
ignacio_b
  • Locked
  • Question
periodic boundary conditions (PBC) using cshift function 1
Replies
1
Views
123
xwb
xwb
peppebig89
  • Locked
  • Question
Hello to all guys, I'm just using t
Replies
4
Views
65
mikrom
mikrom

Part and Inventory Search

Sponsor

Back
Top