jalvesz
commented
9 months ago Hi,
In order to help this idea move forward regarding its extension to sparse algebra I thought about moving parts of FSPARSE here.
I thought about doing a first PR to define the sparse types with something like:
Click to open (stdlib_sparse_kinds.fypp)
```Fortran #:include "common.fypp" !> The `stdlib_sparse_kinds` module provides derived type definitions for different sparse matrices !> !> This code was modified from https://github.com/jalvesz/FSPARSE by its author: Alves Jose module stdlib_sparse_kinds use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp implicit none private ! -- Global parameters enum, bind(C) enumerator :: k_NOSYMMETRY !! Full Sparse matrix (no symmetry considerations) enumerator :: k_SYMTRIINF !! Symmetric Sparse matrix with triangular inferior storage enumerator :: k_SYMTRISUP !! Symmetric Sparse matrix with triangular supperior storage end enum ! -- Classes type, abstract :: sparse_t integer :: nrows = 0 !! number of rows integer :: ncols = 0 !! number of columns integer :: nnz = 0 !! number of non-zero values integer :: sym = k_NOSYMMETRY !! assumed storage symmetry integer :: base = 1 !! index base = 0 for (C) or 1 (Fortran) end type !! COO: COOrdinates compresed format type, public, extends(sparse_t) :: COO_t logical :: isOrdered = .false. !! wether the matrix is ordered or not integer, allocatable :: index(:,:) !! Matrix coordinates index(2,nnz) contains procedure :: malloc => malloc_coo end type #:for k1, t1 in (REAL_KINDS_TYPES) type, public, extends(COO_t) :: COO_${k1}$ ${t1}$, allocatable :: data(:) !! single precision values end type #:endfor !! CSR: Compressed sparse row or Yale format type, extends(sparse_t) :: CSR_t integer, allocatable :: col(:) !! matrix column pointer integer, allocatable :: rowptr(:) !! matrix row pointer contains procedure :: malloc => malloc_csr end type #:for k1, t1 in (REAL_KINDS_TYPES) type, public, extends(CSR_t) :: CSR_${k1}$ ${t1}$, allocatable :: data(:) !! single precision values end type #:endfor !! CSC: Compressed sparse column type, extends(sparse_t) :: CSC_t integer, allocatable :: colptr(:) !! matrix column pointer integer, allocatable :: row(:) !! matrix row pointer contains procedure :: malloc => malloc_csc end type #:for k1, t1 in (REAL_KINDS_TYPES) type, public, extends(CSC_t) :: CSC_${k1}$ ${t1}$, allocatable :: data(:) !! single precision values end type #:endfor !! Compressed ELLPACK type, extends(sparse_t) :: ELL_t integer :: K = 0 !! maximum number of nonzeros per row integer, allocatable :: index(:,:) !! column indices contains procedure :: malloc => malloc_ell end type #:for k1, t1 in (REAL_KINDS_TYPES) type, public, extends(ELL_t) :: ELL_${k1}$ ${t1}$, allocatable :: data(:,:) !! single precision values end type #:endfor contains subroutine malloc_coo(self,num_rows,num_cols,nnz) class(COO_t) :: self integer, intent(in) :: num_rows integer, intent(in) :: num_cols integer, intent(in) :: nnz integer, allocatable :: temp_idx(:,:) !----------------------------------------------------- self%nrows = num_rows self%ncols = num_cols self%nnz = nnz if(.not.allocated(self%index)) then allocate(temp_idx(2,nnz) , source = 0 ) else allocate(temp_idx(2,nnz) , source = self%index ) end if call move_alloc(from=temp_idx,to=self%index) select type(self) #:for k1, t1 in (REAL_KINDS_TYPES) type is(COO_${k1}$) block ${t1}$, allocatable :: temp_data_${k1}$(:) if(.not.allocated(self%data)) then allocate(temp_data_${k1}$(nnz) , source = 0._${k1}$ ) else allocate(temp_data_${k1}$(nnz) , source = self%data ) end if call move_alloc(from=temp_data_${k1}$,to=self%data) end block #:endfor end select end subroutine subroutine malloc_csr(self,num_rows,num_cols,nnz) class(CSR_t) :: self integer, intent(in) :: num_rows integer, intent(in) :: num_cols integer, intent(in) :: nnz integer, allocatable :: temp_idx(:) !----------------------------------------------------- self%nrows = num_rows self%ncols = num_cols self%nnz = nnz if(.not.allocated(self%col)) then allocate(temp_idx(nnz) , source = 0 ) else allocate(temp_idx(nnz) , source = self%col ) end if call move_alloc(from=temp_idx,to=self%col) if(.not.allocated(self%rowptr)) then allocate(temp_idx(num_rows+1) , source = 0 ) else allocate(temp_idx(num_rows+1) , source = self%rowptr ) end if call move_alloc(from=temp_idx,to=self%rowptr) select type(self) #:for k1, t1 in (REAL_KINDS_TYPES) type is(CSR_${k1}$) block ${t1}$, allocatable :: temp_data_${k1}$(:) if(.not.allocated(self%data)) then allocate(temp_data_${k1}$(nnz) , source = 0._${k1}$ ) else allocate(temp_data_${k1}$(nnz) , source = self%data ) end if call move_alloc(from=temp_data_${k1}$,to=self%data) end block #:endfor end select end subroutine subroutine malloc_csc(self,num_rows,num_cols,nnz) class(CSC_t) :: self integer, intent(in) :: num_rows integer, intent(in) :: num_cols integer, intent(in) :: nnz integer, allocatable :: temp_idx(:) !----------------------------------------------------- self%nrows = num_rows self%ncols = num_cols self%nnz = nnz if(.not.allocated(self%row)) then allocate(temp_idx(nnz) , source = 0 ) else allocate(temp_idx(nnz) , source = self%row ) end if call move_alloc(from=temp_idx,to=self%row) if(.not.allocated(self%colptr)) then allocate(temp_idx(num_cols+1) , source = 0 ) else allocate(temp_idx(num_cols+1) , source = self%colptr ) end if call move_alloc(from=temp_idx,to=self%colptr) select type(self) #:for k1, t1 in (REAL_KINDS_TYPES) type is(CSC_${k1}$) block ${t1}$, allocatable :: temp_data_${k1}$(:) if(.not.allocated(self%data)) then allocate(temp_data_${k1}$(nnz) , source = 0._${k1}$ ) else allocate(temp_data_${k1}$(nnz) , source = self%data ) end if call move_alloc(from=temp_data_${k1}$,to=self%data) end block #:endfor end select end subroutine subroutine malloc_ell(self,num_rows,num_cols,num_nz_rows) class(ELL_t) :: self integer, intent(in) :: num_rows !! number of rows integer, intent(in) :: num_cols !! number of columns integer, intent(in) :: num_nz_rows !! number of non zeros per row integer, allocatable :: temp_idx(:,:) !----------------------------------------------------- self%nrows = num_rows self%ncols = num_cols self%K = num_nz_rows if(.not.allocated(self%index)) then allocate(temp_idx(num_rows,num_nz_rows) , source = 0 ) else allocate(temp_idx(num_rows,num_nz_rows) , source = self%index ) end if call move_alloc(from=temp_idx,to=self%index) select type(self) #:for k1, t1 in (REAL_KINDS_TYPES) type is(ELL_${k1}$) block ${t1}$, allocatable :: temp_data_${k1}$(:,:) if(.not.allocated(self%data)) then allocate(temp_data_${k1}$(num_rows,num_nz_rows) , source = 0._${k1}$ ) else allocate(temp_data_${k1}$(num_rows,num_nz_rows) , source = self%data ) end if call move_alloc(from=temp_data_${k1}$,to=self%data) end block #:endfor end select end subroutine end module stdlib_sparse_kinds ```Other formats can be added following a similar pattern. Before going any further with the methods, I though that having some comments about the derived types for the data containment is important.
For instance:
- Should dense matrices be included as well as one of the derived types? or plain allocatables should do?
jvdp1
perazz
gnikit
Hello stdlib developers, I'm opening this issue to summarize and coalesce upcoming efforts to integrate linear algebra operations in stdlib, in particular:
1) Accessible interfaces for common linear algebra operations
We want stdlib to be able to solve linear algebra tasks by wrapping against libraries like BLAS, LAPACK, SCALAPACK, with a user-friendly interface. I think this means:
2) Support for common IO formats for matrices and tensors
To allow efficient working with multidimensional arrays, we should support easy serialization/deserialization and conversion among formats. This could include formats like NPZ, matrixmarket, etc. An easily extendible API should be provided to have plugins for other formats which might or might not fit in the scope of stdlib.
I believe an equally important task is to define derived types that are capable of storing temporary information which is not strictly matrix data (e.g. matrix factorization, or working arrays for the matrix solvers), to avoid unnecessary overhead in case of repeated algebra operations. This means that "simple" array storage may need to be replaced with matrix derived types in those cases.
Because there are plenty of options in defining these APIs, It is crucial to the success of this task that as much feedback as possible is given, so I would like to encourage all ideas - and criticisms - to be discussed on this issue, so that we can come up with the best possible version. I am also opening a discussion page on the Fortran-lang Discourse that we can use for more verbose discussions.
Thank you, Federico
cc @certik @awvwgk @fortran-lang/stdlib @fortran-lang/admins
Linked issues
Linear algebra (BLAS/LAPACK)
1
10
67
450
476
-> Regarding dense algebra, I've started a discussion at https://github.com/fortran-lang/stdlib/issues/450
Sparse algebra
38