stdlib
stdlib copied to clipboard
fortran-lang
Support for linear algebra
Hello stdlib developers, I'm opening this issue to summarize and coalesce upcoming efforts to integrate linear algebra operations in stdlib, in particular:
- 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:
- To get the best possible user-level API (almost no inputs at all), maintaining familiarity the syntax that other libraries also provide (scipy, etc.) but have an expert interface with all control knobs to tune the algorithms (e.g. via settings stored in derived types). Thanks to Fortran arrays, If this is done well, I believe stdlib could have the best linear algebra API out there;
- Wrapping against BLAS and LAPACK means that the low-level API should probably not be changed (so that down the road, stdlib could just link against external BLAS libraries or platform-dependent frameworks), but a reference implementation must be provided anyways. Should we aim at maintaining our own Modernized BLAS/LAPACK at fortran-lang, or just automate the download process from netlib? I would personally like the idea to develop a Modernized version once and forall (those reference implementations will almost never change anymore), although those are huge repos that would require significant work.
- 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
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)
#: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?
@jalvesz it seems a good start to me. I developed a library with similar interfaces. There have been already a lot on discussion about sparse matrices, and it is always difficult to find a clear clonclusion. So, I suggest to keep the first draft a simple as possible. Other formats (and dense matrices) could be added later.
Would it make sense to open a new branch on the stdlib repo to coalesce linear algebra progress around?
Would it make sense to open a new branch on the stdlib repo to coalesce linear algebra progress around?
Yes, it is probably a good idea. There is also #189 related to this topic that is open since a while (cc @jalvesz )
There is also this page related to sparse matrices
Let me know when you open the branch and I'll move the PR under that one, indeed it would be easier to coalesce all linear algebra topics under one branch that could move a bit faster and serve as playground
Hey @perazz is this related to the grant money from STF and the work done over at https://github.com/perazz/fortran-lapack or would this be an independent piece of work (to be performed by the community)? Regardless, it would make sense to have an stdlib branch containing the relevant progress.
