How to resolve error associated with finding the determinant of n x n matrix

Question

I saw this Fortran code to calculate the determinant of a n x n size matrix on the internet. I intend to apply it to my research which involves a very large size of matrix. But I tried to apply it first to a 3 x 3 matrix and run it using gfortran compiler to understand how it works, but I get the following errors:

determinant.f:39.26:

        cf = cofactor(matrix, i, 1)
                      1

Warning: Legacy Extension: REAL array index at (1) determinant.f:39.26:

        cf = cofactor(matrix, i, 1)
                      1

Error: Array index at (1) is an array of rank 2 determinant.f:14.21:

  D = determinant(A)
                 1

Error: Return type mismatch of function 'determinant' at (1) (UNKNOWN/REAL(4)) determinant.f:14.21:

  D = determinant(A)
                 1

Error: Procedure 'determinant' at (1) with assumed-shape dummy argument 'matrix' must have an explicit interface determinant.f:14.10:

  D = determinant(A)
      1

Error: Function 'determinant' at (1) has no IMPLICIT type

The code which I worked on is as shown below:

      Program determinant_matrix

      implicit none
      real, dimension(9) :: Aa
      real, dimension(3,3):: A
      real :: D

      integer, dimension(1:2) :: order2 = (/2,1/)
      Aa = (/2.0,9.0,4.0,7.0,5.0,3.0,6.0,1.0,8.0/)
       A= reshape(Aa,(/3,3/), order = order2 )

      D = determinant(A)

      end program determinant_matrix



   
c       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, cofactor

      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

Please, is there a way to rectify these errors?

Answer 1

It goes without saying that this is ... not a good method. The first row would spawn N cofactors. Each of those would spawn (N-1) cofactors and so on, resulting in N! cofactors. N! (N factorial) grows very fast with N. There are also a lot of allocations/deallocations and there is likely to be a lot of memory in use at the same time.

To use a function taking an assumed-shape array as argument requires an explicit interface; either: a separate module (probably best here), a "contain"-ed internal function, or an explicit interface statement.

You can certainly do the problem by recursion, but you should reduce it first (e.g. by similar row operations to Gauss elimination). Then you end up multiplying the diagonal elements only, and there need be no mention of cofactors.

The following code will do the determinant by recursion, but please note:

(1) it would be far better using an optimised library routine (e.g. Lapack, as suggested by Ian Bush);

(2) it would be far better to use a loop than recursion - there's a real danger of growing memory;

(3) the working will modify matrix A (use a copy if you don't want that);

(4) it's very vulnerable to floating-point round-off error.

module matmodule
   use iso_fortran_env, only: real64
   implicit none
   real(real64), parameter :: TOL = 1.0e-20

   private
   public det

contains

   recursive function det( A ) result( res )
      real(real64) res
      real(real64), intent(inout) :: A(:,:)                ! WARNING! this will change A
      integer N, i, j
      real(real64) temp

      N = size( A, 1 )
      ! Base case for recursion
      if ( N == 1 ) then
         res = A(1,1)
         return
      end if

      res = 1.0
      ! Swap rows (determinant switches sign) so that largest absolute element is at top of first column
      i = 1
      do j = 2, N
         if ( abs( A(j,1) ) > abs( A(i,1) ) ) i = j
      end do
      if ( i /= 1 ) then
         res = -res
         do j = 1, N
            temp = A(1,j);   A(1,j) = A(i,j);   A(i,j) = temp
         end do
      end if

      ! Effectively zero
      if ( abs( A(1,1) ) < TOL ) then
         res = 0.0
         return
      end if

      ! Row operations to eliminate first column
      do i = 2, N
         A(i,:) = A(i,:) - ( A(i,1) / A(1,1) ) * A(1,:)
      end do

      ! Expand determinant on lone element in column 1
      res = res * A(1,1) * det( A(2:N,2:N) )

   end function det

end module matmodule

!=======================================================================

program test
   use iso_fortran_env, only: real64
   use matmodule, only: det
   implicit none
   real(real64) :: A(3,3) = reshape( [ 2.0, 9.0, 4.0, 7.0, 5.0, 3.0, 6.0, 1.0, 8.0 ], [ 3, 3 ], order = [ 2, 1 ] )

   print *, det( A )

end program test

Output:

  -360.00000000000000