Gauss Algorithm written in Fortran slower than written in Python for large matrices?

Help
9

For row-intensive operations, like Gauss elimination etc., one could possibly prepare a hybrid matrix, constructed as an array of pointers (embedded in a derived type, as Fortran, lamentably, does not provide regular array of pointers) to the matrix rows. That way, the row data are contiguous and all operations are done on rank-1 arrays.
Due to the aforementioned deficiency, the syntax is a bit clumsy, but it works, one can swap whole or partial rows as seen in the following super simple sample.

program main
  implicit none
  type rows
    real, pointer :: row(:)
  end type rows
  type(rows), allocatable :: array(:)
  integer :: n=8, i, j
  real, allocatable :: swap(:)
  real :: x,y

  allocate(array(n),swap(n))
  do i=1,n
    allocate(array(i)%row(n))
    do j=1,n
      array(i)%row(j) = 100.0*i+j
    end do
  end do
  print 10,array(1)%row  !  101.0 102.0 103.0 104.0 105.0 106.0 107.0 108.0
  print 10,array(5)%row  !  501.0 502.0 503.0 504.0 505.0 506.0 507.0 508.0
  swap = array(1)%row
  array(1)%row = array(5)%row
  array(5)%row = swap
  print 10,array(1)%row  !  501.0 502.0 503.0 504.0 505.0 506.0 507.0 508.0
  print 10,array(5)%row  !  101.0 102.0 103.0 104.0 105.0 106.0 107.0 108.0
  swap(1:n/2) = array(1)%row(1:n/2)
  array(1)%row(1:n/2) = array(5)%row(1:n/2)
  array(5)%row(1:n/2) = swap(1:n/2)
  print 10,array(1)%row  !  101.0 102.0 103.0 104.0 505.0 506.0 507.0 508.0
  print 10,array(5)%row  !  501.0 502.0 503.0 504.0 105.0 106.0 107.0 108.0
10 format(*(f6.1))
end program main

BTW, anybody knowing why actually array of pointers cannot exist in Modern Fortran?

1 Like