Comments (6)
jalvesz
commented on June 12, 2024
1
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_kindsOther 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?
from stdlib.
jvdp1
commented on June 12, 2024
1
@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.
from stdlib.
gnikit
commented on June 12, 2024
1
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.
from stdlib.
perazz
commented on June 12, 2024
Would it make sense to open a new branch on the stdlib repo to coalesce linear algebra progress around?
from stdlib.
jvdp1
commented on June 12, 2024
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
from stdlib.
jalvesz
commented on June 12, 2024
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
from stdlib.
Related Issues (20)
- hash_functions test fails on i386: `Segmentation fault - invalid memory reference`
- Request to upgrade Intel-classic compiler in macOS CI
- Add `library` configuration to `stdlib-fpm`
- Massive slow down in docs generation HOT 4
- Unexpected performance of hash maps HOT 8
- python preprocessor HOT 11
- add topic tags `lapack`, `blas`, `linear-algebra` HOT 1
- Improve descriptions of rotm, rotmg, stdlib_srotm, stdlib_srotmg
- Don't repeat names of procedures in descriptions
- stdlib_io_npy, FPM and Rank > 4 HOT 5
- Missing documentation for `LAPACK`-related functions HOT 6
- Adding the logarithmic derivative of the gamma function (digamma) to stdlib_specialfunctions_gamma HOT 1
- example_starts_with prints logical results as binary
- Nonstandard forward reference to 'lk' is not allowed in the same specification part causes compilation errors HOT 4
- Issue/Question about the output of `lstsq` HOT 3
- Building and testing stdlib on WSL1 with two fails HOT 5
- Adding matrix norms to the `stdlib_linalg` module. HOT 8
- Extend `sort_index` interface to allow `int32` index argument HOT 1
- Bug in the complex least-squares solver HOT 6
- Conflict between the pure function `stdlib_iparam2stage` and `openmp` HOT 1
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from stdlib.