Skip to content

Stdlib Sparse matrix API

Jump to bottom
Williams A. Lima edited this page May 20, 2020 · 43 revisions

Stdlib Sparse matrix API Specification

https://github.com/fortran-lang/stdlib/issues/38

Introduction

Goals

To be compact instead of being exhaustive. It aims at supplying Fortran users with a minimum (yet useful) number of routines and data structures related to sparse matrices storage and operations. This library is particularly targeted at a non-expert in numerical computation public. Thus we aim at having a simple and easy to use API.

Sparse matrix representations supported

This section is based on Saad (1994). In that work, a much more complete and extensive list of formats is described. Here we take only the ones that we think are most useful at the moment.

Coordinate format (COO)

Given an m by n real or complex matrix A containing nnz nonzero elements with each element denoted by a_ij this format represents A using a set of three arrays: values, rows, and columns, as described below.

values A real/complex array of size nnz containing the matrix elements a_ij in any order.

rows An integer array of size nnz containing the row indices of the elements a_ij.

columns An integer array of size nnz containing the column indices of the elements a_ij.

The Sparse_COO_t type

type :: Sparse_Matrix_COO_t
   real (kind = prec), allocatable, dimension (:) :: values
   integer, allocatable, dimension (:) :: rows
   integer, allocatable, dimension (:) :: columns
end type Sparse_Matrix_COO_t

Attributes

**Need to discuss sizes for integer.

values real, dynamically allocated, rank 1 array.

The non-zero values of the matrix A in any order.

rows Integer, dynamically allocated, rank 1 array.

The corresponding row indices of the values stored in the array values.

columns Integer, dynamically allocated, rank 1 array.

The corresponding column indices of the values stored in the arra values.

Compressed Sparse Row (CSR)

Given an m by n real or complex matrix A containing nnz nonzero elements with each element denoted by a_ij this format represents A using a set of three arrays: values, ja, and ia, as described below.

values A real/complex array of size nnz containing the matrix elements a_ij stored row by row from row 1 to row n.

ja An integer array of size nnz containing the column indices of the elements a_ij as stored in the array values.

ia An integer array of size n+1 containing the index in the arrays values and ja where each row starts. The value at ia(n+1) always has the value ia(1)+nnz.

The Sparse_CSR_t type

type :: Sparse_Matrix_CSR_t
   real (kind = prec), allocatable, dimension (:) :: values
   integer, allocatable, dimension (:) :: ja
   integer, allocatable, dimension (:) :: ia
end type Sparse_Matrix_COO_t

Attributes

**Need to discuss sizes for integer.

values real, dynamically allocated, rank 1 array.

The non-zero values of the matrix A in any order.

ja Integer, dynamically allocated, rank 1 array.

The column indices of the elements a_ij.

ia Integer, dynamically allocated, rank 1 array.

The indices in the arrays values and ja where each row starts.

Compressed Sparse Column (CSC)

This format is similar to the CSR format described previously. The difference is that instead of storing row values we store column values in the array values. The exact description of this format is given below.

Given an m by n real or complex matrix A containing nnz nonzero elements with each element denoted by a_ij this format represents A using a set of three arrays: values, ja, and ia, as described below.

values A real/complex array of size nnz containing the matrix elements a_ij stored column by column from column 1 to column m.

ia An integer array of size nnz containing the row indices of the elements a_ij as stored in the array values.

ja An integer array of size m+1 containing the index in the arrays values and ia where each column starts. The value at ja(m+1) always has the value ja(1)+nnz.

The Sparse_CSC_t type

type :: Sparse_Matrix_CSR_t
   real (kind = prec), allocatable, dimension (:) :: values
   integer, allocatable, dimension (:) :: ja
   integer, allocatable, dimension (:) :: ia
end type Sparse_Matrix_COO_t

Attributes

**Need to discuss sizes for integer.

values real, dynamically allocated, rank 1 array.

The non-zero values of the matrix A in any order.

ja Integer, dynamically allocated, rank 1 array.

The indices in the arrays values and ia where each column starts.

ia Integer, dynamically allocated, rank 1 array.

The row indices of the elements a_ij.

Creational subroutines

Sparse_matrix_create_coo

The subroutine Sparse_matrix_create_coo initializes a variable of type Sparse_Matrix_COO_t.

Syntax

CALL Sparse_matrix_create_coo (A, m, n, nnz, err)

Arguments

    A         The new matrix to be created.

    m, n     Integers, scalars corresponding to the dimensions of A.

    nnz      Integer, scalar, the number of non-zero entries in A.

    err        Integer, scalar, 0 if no error occurred.

Example

program Test_sparse
   type (Sparse_Matrix_COO_t) :: A
   integer :: ierr
   integer :: m, n, nnz

   m = 1_090_920
   n = 1_090_920
   nnz = 3_083_796

   call Sparse_matrix_create_coo (A, m, n, nnz, ierr)

end program Test_sparse

Conversion subroutines

Algebraic operations

All the core routines work only with the COO format. All their input/output arguments are in the COO format. On the other hand, the use of the supplied conversion routines and a generic interface makes the use of the core routines transparent for the end-user. fig01

Level 1 - Vector-vector operations

Function Name

Syntax

Arguments

Return value

Description

Level 2 - Matrix-vector operations

Level 3 - Matrix-matrix operations

Sparse_matmul_coo (Low level)

The subroutine Sparse_matmul_coo performs matrix multiplication of two COO sparse matrices.

The input values A_i, A_j, A_values and B_i, B_j, B_values must be previously initialized and have the correct dimensions. The output matrix variables (C_i, C_j, C_values) are initialized inside Sparse_matmul_coo, if already initialized when calling this subroutine, then on output all previous information in them will have been overwritten.

Syntax

CALL Sparse_matmul_coo (A_i, A_j, A_values, B_i, B_j, B_values, C_i, C_j, C_values, err)

Arguments

    A_i     Row indices of input matrix A.

    A_j     Column indices of input matrix A.

    A_values     Nonzero values of input matrix A.

    B_i     Row indices of input matrix B.

    B_j     Column indices of input matrix B.

    B_values     Nonzero values of input matrix B.

    C_i     Row indices of input matrix C.

    C_j     Column indices of input matrix C.

    C_values     Nonzero values of input matrix C.

    err    Integer, scalar, 0 if no error occurred.

Example

program Test_sparse
   integer, allocatable, dimension (:) :: A_i, A_j, B_i, B_j, C_i, C_j
   real (kind = prec), allocatable, dimension (:) :: A_values, B_values, C_values
   integer :: ierr
   integer :: m, n, nnz

   m = 100; n = 300; nnz = 952
   allocate (A_i(nnz), A_j(nnz), A_values(nnz))

   m = 300; n = 500; nnz = 1320
   allocate (B_i(nnz), B_j(nnz), B_values(nnz))

   !...Populate matrices A and B

   call Sparse_matmul_coo (A_i, A_j, A_values, B_i, B_j, B_values, C_i, C_j, C_values, ierr) ! Matrix C is initialized in Sparse_matmul_coo

end program Test_sparse

Sparse_matmul_coo (Middle level)

The subroutine Sparse_matmul_coo performs matrix multiplication of two Sparse_Matrix_COO_t.

The input matrices A and B must be previously initialized and have the correct dimensions. The output matrix C is initialized inside Sparse_matmul_coo. If C is already initialized when calling this subroutine, then on output all previous information in it will have been overwritten.

Syntax

CALL Sparse_matmul_coo (A, B, C, err)

Arguments

    A     A Sparse_Matrix_COO_t.

    B     A Sparse_Matrix_COO_t.

    C     A Sparse_Matrix_COO_t.

    err    Integer, scalar, 0 if no error occurred.

Implementation

subroutine Sparse_matmul_coo (A, B, C)
   type (Sparse_Matrix_COO_t), intent (in)  :: A, B
   type (Sparse_Matrix_COO_t), intent (out) :: C

   ! Call low level Sparse_matmul_coo
   ! Matrix C is initialized in Sparse_matmul_coo
   call Sparse_matmul_coo (A%rows, A%columns, A%values, &
                           B%rows, B%columns, B%values, &
                           C%rows, C%columns, C%values, ierr)
end subroutine Sparse_matmul_coo

Example

program Test_sparse
   type (Sparse_Matrix_COO_t) :: A, B, C
   integer :: ierr
   integer :: m, n, nnz

   m = 100
   n = 300
   nnz = 952
   call Sparse_matrix_create_coo (A, m, n, nnz, ierr)

   m = 300
   n = 500
   nnz = 1320
   call Sparse_matrix_create_coo (A, m, n, nnz, ierr)

   !...Populate matrices A and B

   ! Matrix C is initialized in Sparse_matmul_coo
   call Sparse_matmul_coo (A, B, C, ierr)

end program Test_sparse

Sparse_matmul_csr (Low level)

The subroutine Sparse_matmul_csr performs matrix multiplication of two CSR matrices (A, B). Actually this subroutine is written using Sparse_matmul_coo. On input, matrices A and B are converted to COO representation, then Sparse_matmul_coo is called, then the result is converted back to CSR format.

The input matrices A and B must be previously initialized and have the correct dimensions. The output matrix C is initialized inside Sparse_matmul_csr. If C is already initialized when calling this subroutine, then on output all previous information in it will have been overwritten.

Syntax

CALL Sparse_matmul_csr (A_ja, A_ia, A_values, B_ja, B_ia, B_values, C_ja, C_ia, C_values, err)

Arguments

    A_ja     Column indices of the elements a_ij.

    A_ia     Index in the arrays values and A_ja where each row starts.

    A_values     Contains the matrix elements a_ij stored row by row from row 1 to row n.

    B_ja     Column indices of the elements b_ij.

    B_ia     Index in the arrays values and B_ja where each row starts.

    B_values     Contains the matrix elements b_ij stored row by row from row 1 to row n.

    err    Integer, scalar, 0 if no error occurred.

    C_ia     Index in the arrays values and C_ja where each row starts.

    C_values     Contains the matrix elements c_ij stored row by row from row 1 to row n.

    err    Integer, scalar, 0 if no error occurred.

Implementation

subroutine sparse_matmul_csr (A_ja, A_ia, A_values, B_ja, B_ia, B_values, C_ja, C_ia, C_values, ierr)
  ! Arguments
  integer, dimension (:), intent (in)  :: A_ja, A_ia, B_ja, B_ia
  integer, allocatable, dimension (:), intent (out) :: C_ja, C_ia
  real (kind = prec), dimension (:), intent (in)  :: A_values, B_values
  real (kind = prec), allocatable, dimension (:), intent (out) :: C_values
  ! Local variables
  integer, allocatable, dimension (:) :: A_rows_coo, A_cols_coo
  integer, allocatable, dimension (:) :: B_rows_coo, B_cols_coo
  integer, allocatable, dimension (:) :: C_rows_coo, C_cols_coo
  real (kind = prec), allocatable, dimension (:) :: A_values_coo, B_values_coo, C_values_coo

  call Sparse_Matrix_CSR_to_COO (A_ja, A_ia, A_values, A_rows_coo, A_cols_coo, A_values_coo)
  call Sparse_Matrix_CSR_to_COO (B_ja, B_ia, B_values, B_rows_coo, B_cols_coo, B_values_coo)

  call Sparse_matmul (A_rows_coo, A_cols_coo, A_values_coo, &
                      B_rows_coo, B_cols_coo, B_values_coo, &
                      C_rows_coo, C_cols_coo, C_values_coo, ierr)
  
  call Sparse_Matrix_COO_to_CSR (C_rows_coo, C_cols_coo, C_values_coo, C_ja, C_ia, C_values, ierr)

end subroutine sparse_matmul_csr

Sparse_matmul_csr (Middle level)

The subroutine Sparse_matmul_csr performs matrix multiplication of two Sparse_Matrix_CSR_t. Actually this subroutine is written using Sparse_matmul_coo. On input, matrices A and B are converted to Sparse_matmul_coo_t, then Sparse_matmul_coo es called and the result is converted back to Sparse_Matrix_CSR_t.

The input matrices A and B must be previously initialized and have the correct dimensions. The output matrix C is initialized inside Sparse_matmul_csr. If C is already initialized when calling this subroutine, then on output all previous information in it will have been overwritten.

Syntax

CALL Sparse_matmul_csr (A, B, C, err)

Arguments

    A     A Sparse_Matrix_CSR_t.

    B     A Sparse_Matrix_CSR_t.

    C     A Sparse_Matrix_CSR_t.

    err    Integer, scalar, 0 if no error occurred.

Implementation

subroutine sparse_matmul_csr (a, b, c)
  ! Arguments
  type (Sparse_Matrix_CSR_t), intent (in)  :: a, b
  type (Sparse_Matrix_CSR_t), intent (out) :: c
  ! Local variables
  type (Sparse_Matrix_COO_t) :: a_coo, b_coo, c_coo

  call Sparse_Matrix_CSR_to_COO (a, a_coo)
  call Sparse_Matrix_CSR_to_COO (b, b_coo)

  call Sparse_matmul (a_coo, b_coo, c_coo)
  
  call Sparse_Matrix_COO_to_CSR (c_coo, c)

end subroutine sparse_matmul_csr

Sparse_matmul_csc

The subroutine Sparse_matmul_csc performs matrix multiplication of two Sparse_Matrix_CSC_t. Actually this subroutine is written using Sparse_matmul_coo. On input, matrices A and B are converted to Sparse_matmul_coo_t, then Sparse_matmul_coo es called and the result is converted back to Sparse_Matrix_CSC_t.

The input matrices A and B must be previously initialized and have the correct dimensions. The output matrix C is initialized inside Sparse_matmul_csr. If C is already initialized when calling this subroutine, then on output all previous information in it will have been overwritten.

Syntax

CALL Sparse_matmul_csc (A, B, C, err)

Arguments

    A     A Sparse_Matrix_CSC_t.

    B     A Sparse_Matrix_CSC_t.

    C     A Sparse_Matrix_CSC_t.

    err    Integer, scalar, 0 if no error occurred.

Implementation

subroutine sparse_matmul_csc (a, b, c)
  ! Arguments
  type (Sparse_Matrix_CSC_t), intent (in)  :: a, b
  type (Sparse_Matrix_CSC_t), intent (out) :: c
  ! Local variables
  type (Sparse_Matrix_COO_t) :: a_coo, b_coo, c_coo

  call Sparse_Matrix_CSR_to_CSC (a, a_coo)
  call Sparse_Matrix_CSR_to_CSC (b, b_coo)

  call Sparse_matmul (a_coo, b_coo, c_coo)
  
  call Sparse_Matrix_COO_to_CSC (c_coo, c)

end subroutine sparse_matmul_csc

Generic interface for Sparse_matmul

interface sparse_matmul
   module procedure sparse_matmul_csr
   module procedure sparse_matmul_ccr
   module procedure sparse_matmul_coo
end interface sparse_matmul

Utilities

Input/Output

Sorting

References

Saad, Y., SPARSKIT: A basic tool kit for sparse matrix computation, 1994.https://www-users.cs.umn.edu/~saad/software/SPARSKIT/

Clone this wiki locally