Sparse array in Fortran

Question

I am very new to Fortran (a month or so). Before I have been mostly scripting in Python. I am currently working on a code in which I need some method to store a sparse array. The array is square and would be roughly 3E6-by-3E6. At first I tried to store the complete array with 8-bit reals but encountered "insufficient virtual memory".

To reduce memory consumption I wrote a set of very simple routines and a type to store the matrix as a row, column and value triple. By using another array to copy to and then deallocate and allocate I could store values this triple format. The code turned out dead slow. Of course it slowed down as a function of how many values I stored in the sparse array.

I noticed that there is a very simple way to reallocate in Fortran 2003 and onwards using vec = [vec(:), val] to append one element. This requires compilation with -assume realloc_lhs. Unfortunately this is exactly as slow as the version I wrote before.

One problem is that although I know how large the matrix would be, I only know the upper limit on the number of non-zero elements per row (or column). I haven't exploited that information and that is why I try to reallocate.

I have attached the current versions that uses -assume realloc_lhs. I am grateful for any kind of tip or hint towards a better solution.

   TYPE sprs
      INTEGER               :: n, len
      REAL,    ALLOCATABLE  :: val(:)
      INTEGER, ALLOCATABLE  :: irow(:)
      INTEGER, ALLOCATABLE  :: icol(:)
   END TYPE sprs

!===================================================================================================
   SUBROUTINE store_sprs(val,irow,icol,matrix)
!===================================================================================================
!
   IMPLICIT NONE
!
! Arguments:
!-----------

   REAL,    INTENT(IN)       :: val
   INTEGER, INTENT(IN)       :: irow
   INTEGER, INTENT(IN)       :: icol
   TYPE(sprs), INTENT(INOUT) :: matrix
!
! Locals:
!--------

   IF (ABS(val)>=1.0E-50) THEN
      CALL add2list_re(val, matrix.val)
      CALL add2list_int(irow, matrix.irow)
      CALL add2list_int(icol, matrix.icol)
   ENDIF
   END SUBROUTINE store_sprs


   SUBROUTINE add2list_re(val,vec)
!===================================================================================================
   IMPLICIT NONE
!
! Arguments:
!-----------

   REAL, INTENT(IN) :: val
   REAL, ALLOCATABLE, INTENT(INOUT) :: vec(:)

!
! Locals:
!--------

   vec = [vec(:), val]

   END SUBROUTINE add2list_re

!===================================================================================================
   SUBROUTINE add2list_int(val,vec)
!===================================================================================================
   IMPLICIT NONE
!
! Arguments:
!-----------

   INTEGER, INTENT(IN) :: val
   INTEGER, ALLOCATABLE, INTENT(INOUT) :: vec(:)

!
! Locals:
!--------

   vec = [vec(:), val]

   END SUBROUTINE add2list_int

Answer 1

There are a number of data structures which handle sparse data well. I have found the C++ std::Vector model to work nicely.

This model allocates space for more elements than are needed and keeps track of how many elements are actually in use. Most of the time when an element is added, no new memory allocation is needed. When more space is needed, the allocated space is doubled. This leads to O(log(n)) allocations and much better performance.

This can be implemented in Fortran using a matrix class and a matrix element class, for example,

module sparse_matrix
  implicit none

  private
  public :: dp
  public :: SparseElement
  public :: Sparse

  integer, parameter :: dp=selected_real_kind(15,300)

  type SparseElement
    integer  :: irow
    integer  :: icol
    real(dp) :: val
  end type

  type Sparse
    ! Only the first no_elements elements will be in use.
    integer                          :: no_elements
    type(SparseElement), allocatable :: elements_(:)
  contains
    procedure, public :: elements
    procedure, public :: add_element
  end type

  interface Sparse
    module procedure new_Sparse
  end interface
contains

! Initialise the Sparse matrix to an empty matrix.
function new_Sparse() result(this)
  implicit none

  type(Sparse) :: this

  this%no_elements = 0
  allocate(this%elements_(0))
end function

! Return only the elements which are in use.
function elements(this) result(output)
  implicit none

  class(Sparse), intent(in)        :: this
  type(SparseElement), allocatable :: output(:)

  output = this%elements_(:this%no_elements)
end function

! Add an element to the array, automatically allocating memory if needed.
subroutine add_element(this,element)
  implicit none

  class(Sparse),       intent(inout) :: this
  type(SparseElement), intent(in)    :: element

  type(SparseElement), allocatable :: temp(:)

  this%no_elements = this%no_elements+1

  ! This is where memory is added.
  ! If this%elements_ would overflow, its length is doubled.
  if (this%no_elements>size(this%elements_)) then
    allocate(temp(2*this%no_elements))
    temp(:this%no_elements-1) = this%elements_
    this%elements_ = temp
  endif

  this%elements_(this%no_elements) = element
end subroutine
end module

program main
  use sparse_matrix
  implicit none

  type(Sparse) :: matrix

  type(SparseElement), allocatable :: elements(:)

  integer :: i

  ! Initialise the matrix.
  matrix = Sparse()

  ! Add some elements to the matrix.
  call matrix%add_element(SparseElement(1,1,1.0_dp))
  call matrix%add_element(SparseElement(3,5,7.0_dp))
  call matrix%add_element(SparseElement(100,200,300.0_dp))

  ! Retrieve the elements from the matrix.
  elements = matrix%elements()

  do i=1,size(elements)
    write(*,*) elements(i)%irow, elements(i)%icol, elements(i)%val
  enddo
end program