Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 3369013

Browse files
authoredNov 9, 2024··
Sparse algebra support with OOP API (#760)
2 parents 427bc68 + 22cd23e commit 3369013

14 files changed

+3079
-0
lines changed
 

‎doc/specs/stdlib_sparse.md

Lines changed: 364 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,364 @@
1+
---
2+
title: sparse
3+
---
4+
5+
# The `stdlib_sparse` module
6+
7+
[TOC]
8+
9+
## Introduction
10+
11+
The `stdlib_sparse` module provides derived types for standard sparse matrix data structures. It also provides math kernels such as sparse matrix-vector product and conversion between matrix types.
12+
13+
## Sparse matrix derived types
14+
15+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
16+
### The `sparse_type` abstract derived type
17+
#### Status
18+
19+
Experimental
20+
21+
#### Description
22+
The parent `sparse_type` is as an abstract derived type holding the basic common meta data needed to define a sparse matrix, as well as shared APIs. All sparse matrix flavors are extended from the `sparse_type`.
23+
24+
```Fortran
25+
type, public, abstract :: sparse_type
26+
integer :: nrows !! number of rows
27+
integer :: ncols !! number of columns
28+
integer :: nnz !! number of non-zero values
29+
integer :: storage !! assumed storage symmetry
30+
end type
31+
```
32+
33+
The storage integer label should be assigned from the module's internal enumerator containing the following three enums:
34+
35+
```Fortran
36+
enum, bind(C)
37+
enumerator :: sparse_full !! Full Sparse matrix (no symmetry considerations)
38+
enumerator :: sparse_lower !! Symmetric Sparse matrix with triangular inferior storage
39+
enumerator :: sparse_upper !! Symmetric Sparse matrix with triangular supperior storage
40+
end enum
41+
```
42+
In the following, all sparse kinds will be presented in two main flavors: a data-less type `<matrix>_type` useful for topological graph operations. And real/complex valued types `<matrix>_<kind>_type` containing the `data` buffer for the matrix values. The following rectangular matrix will be used to showcase how each sparse matrix holds the data internally:
43+
44+
$$ M = \begin{bmatrix}
45+
9 & 0 & 0 & 0 & -3 \\
46+
4 & 7 & 0 & 0 & 0 \\
47+
0 & 8 & -1 & 8 & 0 \\
48+
4 & 0 & 5 & 6 & 0 \\
49+
\end{bmatrix} $$
50+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
51+
### `COO`: The COOrdinates compressed sparse format
52+
#### Status
53+
54+
Experimental
55+
56+
#### Description
57+
The `COO`, triplet or `ijv` format defines all non-zero elements of the matrix by explicitly allocating the `i,j` index and the value of the matrix. While some implementations use separate `row` and `col` arrays for the index, here we use a 2D array in order to promote fast memory acces to `ij`.
58+
59+
```Fortran
60+
type(COO_sp_type) :: COO
61+
call COO%malloc(4,5,10)
62+
COO%data(:) = real([9,-3,4,7,8,-1,8,4,5,6])
63+
COO%index(1:2,1) = [1,1]
64+
COO%index(1:2,2) = [1,5]
65+
COO%index(1:2,3) = [2,1]
66+
COO%index(1:2,4) = [2,2]
67+
COO%index(1:2,5) = [3,2]
68+
COO%index(1:2,6) = [3,3]
69+
COO%index(1:2,7) = [3,4]
70+
COO%index(1:2,8) = [4,1]
71+
COO%index(1:2,9) = [4,3]
72+
COO%index(1:2,10) = [4,4]
73+
```
74+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
75+
### `CSR`: The Compressed Sparse Row or Yale format
76+
#### Status
77+
78+
Experimental
79+
80+
#### Description
81+
The Compressed Sparse Row or Yale format `CSR` stores the matrix structure by compressing the row indices with a counter pointer `rowptr` enabling to know the first and last non-zero column index `col` of the given row.
82+
83+
```Fortran
84+
type(CSR_sp_type) :: CSR
85+
call CSR%malloc(4,5,10)
86+
CSR%data(:) = real([9,-3,4,7,8,-1,8,4,5,6])
87+
CSR%col(:) = [1,5,1,2,2,3,4,1,3,4]
88+
CSR%rowptr(:) = [1,3,5,8,11]
89+
```
90+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
91+
### `CSC`: The Compressed Sparse Column format
92+
#### Status
93+
94+
Experimental
95+
96+
#### Description
97+
The Compressed Sparse Colum `CSC` is similar to the `CSR` format but values are accesed first by column, thus an index counter is given by `colptr` which enables to know the first and last non-zero row index of a given colum.
98+
99+
```Fortran
100+
type(CSC_sp_type) :: CSC
101+
call CSC%malloc(4,5,10)
102+
CSC%data(:) = real([9,4,4,7,8,-1,5,8,6,-3])
103+
CSC%row(:) = [1,2,4,2,3,3,4,3,4,1]
104+
CSC%colptr(:) = [1,4,6,8,10,11]
105+
```
106+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
107+
### `ELLPACK`: ELL-pack storage format
108+
#### Status
109+
110+
Experimental
111+
112+
#### Description
113+
The `ELL` format stores data in a dense matrix of $nrows \times K$ in column major order. By imposing a constant number of elements per row $K$, this format will incur in additional zeros being stored, but it enables efficient vectorization as memory acces is carried out by constant sized strides.
114+
115+
```Fortran
116+
type(ELL_sp_type) :: ELL
117+
call ELL%malloc(num_rows=4,num_cols=5,num_nz_row=3)
118+
ELL%data(1,1:3) = real([9,-3,0])
119+
ELL%data(2,1:3) = real([4,7,0])
120+
ELL%data(3,1:3) = real([8,-1,8])
121+
ELL%data(4,1:3) = real([4,5,6])
122+
123+
ELL%index(1,1:3) = [1,5,0]
124+
ELL%index(2,1:3) = [1,2,0]
125+
ELL%index(3,1:3) = [2,3,4]
126+
ELL%index(4,1:3) = [1,3,4]
127+
```
128+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
129+
### `SELL-C`: The Sliced ELLPACK with Constant blocks format
130+
#### Status
131+
132+
Experimental
133+
134+
#### Description
135+
The Sliced ELLPACK format `SELLC` is a variation of the `ELLPACK` format. This modification reduces the storage size compared to the `ELLPACK` format but maintaining its efficient data access scheme. It can be seen as an intermediate format between `CSR` and `ELLPACK`. For more details read [the reference](https://arxiv.org/pdf/1307.6209v1)
136+
137+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
138+
## `add`- sparse matrix data accessors
139+
140+
### Status
141+
142+
Experimental
143+
144+
### Description
145+
Type-bound procedures to enable adding data in a sparse matrix.
146+
147+
### Syntax
148+
149+
`call matrix%add(i,j,v)` or
150+
`call matrix%add(i(:),j(:),v(:,:))`
151+
152+
### Arguments
153+
154+
`i`: Shall be an integer value or rank-1 array. It is an `intent(in)` argument.
155+
156+
`j`: Shall be an integer value or rank-1 array. It is an `intent(in)` argument.
157+
158+
`v`: Shall be a `real` or `complex` value or rank-2 array. The type shall be in accordance to the declared sparse matrix object. It is an `intent(in)` argument.
159+
160+
## `at`- sparse matrix data accessors
161+
162+
### Status
163+
164+
Experimental
165+
166+
### Description
167+
Type-bound procedures to enable requesting data from a sparse matrix.
168+
169+
### Syntax
170+
171+
`v = matrix%at(i,j)`
172+
173+
### Arguments
174+
175+
`i` : Shall be an integer value. It is an `intent(in)` argument.
176+
177+
`j` : Shall be an integer value. It is an `intent(in)` argument.
178+
179+
`v` : Shall be a `real` or `complex` value in accordance to the declared sparse matrix object. If the `ij` tuple is within the sparse pattern, `v` contains the value in the data buffer. If the `ij` tuple is outside the sparse pattern, `v` is equal `0`. If the `ij` tuple is outside the matrix pattern `(nrows,ncols)`, `v` is `NaN`.
180+
181+
## Example
182+
```fortran
183+
{!example/linalg/example_sparse_data_accessors.f90!}
184+
```
185+
186+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
187+
## `spmv` - Sparse Matrix-Vector product
188+
189+
### Status
190+
191+
Experimental
192+
193+
### Description
194+
195+
Provide sparse matrix-vector product kernels for the current supported sparse matrix types.
196+
197+
$$y=\alpha*op(M)*x+\beta*y$$
198+
199+
### Syntax
200+
201+
`call ` [[stdlib_sparse_spmv(module):spmv(interface)]] `(matrix,vec_x,vec_y [,alpha,beta,op])`
202+
203+
### Arguments
204+
205+
`matrix`: Shall be a `real` or `complex` sparse type matrix. It is an `intent(in)` argument.
206+
207+
`vec_x`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. It is an `intent(in)` argument.
208+
209+
`vec_y`: Shall be a rank-1 or rank-2 array of `real` or `complex` type array. . It is an `intent(inout)` argument.
210+
211+
`alpha`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `alpha=1`. It is an `intent(in)` argument.
212+
213+
`beta`, `optional` : Shall be a scalar value of the same type as `vec_x`. Default value `beta=0`. It is an `intent(in)` argument.
214+
215+
`op`, `optional`: In-place operator identifier. Shall be a `character(1)` argument. It can have any of the following values: `N`: no transpose, `T`: transpose, `H`: hermitian or complex transpose. These values are provided as constants by the `stdlib_sparse` module: `sparse_op_none`, `sparse_op_transpose`, `sparse_op_hermitian`
216+
217+
<!-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -->
218+
## Sparse matrix to matrix conversions
219+
220+
### Status
221+
222+
Experimental
223+
224+
### Description
225+
226+
This module provides facility functions for converting between storage formats.
227+
228+
### Syntax
229+
230+
`call ` [[stdlib_sparse_conversion(module):coo2ordered(interface)]] `(coo[,sort_data])`
231+
232+
### Arguments
233+
234+
`COO` : Shall be any `COO` type. The same object will be returned with the arrays reallocated to the correct size after removing duplicates. It is an `intent(inout)` argument.
235+
236+
`sort_data`, `optional` : Shall be a `logical` argument to determine whether data in the COO graph should be sorted while sorting the index array, default `.false.`. It is an `intent(in)` argument.
237+
238+
### Syntax
239+
240+
`call ` [[stdlib_sparse_conversion(module):from_ijv(interface)]] `(sparse,row,col[,data,nrows,ncols,num_nz_rows,chunk])`
241+
242+
### Arguments
243+
244+
`sparse` : Shall be a `COO`, `CSR`, `ELL` or `SELLC` type. The graph object will be returned with a canonical shape after sorting and removing duplicates from the `(row,col,data)` triplet. If the graph is `COO_type` no data buffer is allowed. It is an `intent(inout)` argument.
245+
246+
`row` : rows index array. It is an `intent(in)` argument.
247+
248+
`col` : columns index array. It is an `intent(in)` argument.
249+
250+
`data`, `optional`: `real` or `complex` data array. It is an `intent(in)` argument.
251+
252+
`nrows`, `optional`: number of rows, if not given it will be computed from the `row` array. It is an `intent(in)` argument.
253+
254+
`ncols`, `optional`: number of columns, if not given it will be computed from the `col` array. It is an `intent(in)` argument.
255+
256+
`num_nz_rows`, `optional`: number of non zeros per row, only valid in the case of an `ELL` matrix, by default it will computed from the largest row. It is an `intent(in)` argument.
257+
258+
`chunk`, `optional`: chunk size, only valid in the case of a `SELLC` matrix, by default it will be taken from the `SELLC` default attribute chunk size. It is an `intent(in)` argument.
259+
260+
## Example
261+
```fortran
262+
{!example/linalg/example_sparse_from_ijv.f90!}
263+
```
264+
### Syntax
265+
266+
`call ` [[stdlib_sparse_conversion(module):diag(interface)]] `(matrix,diagonal)`
267+
268+
### Arguments
269+
270+
`matrix` : Shall be a `dense`, `COO`, `CSR` or `ELL` type. It is an `intent(in)` argument.
271+
272+
`diagonal` : A rank-1 array of the same type as the `matrix`. It is an `intent(inout)` and `allocatable` argument.
273+
274+
#### Note
275+
If the `diagonal` array has not been previously allocated, the `diag` subroutine will allocate it using the `nrows` of the `matrix`.
276+
277+
### Syntax
278+
279+
`call ` [[stdlib_sparse_conversion(module):dense2coo(interface)]] `(dense,coo)`
280+
281+
### Arguments
282+
283+
`dense` : Shall be a rank-2 array of `real` or `complex` type. It is an `intent(in)` argument.
284+
285+
`coo` : Shall be a `COO` type of `real` or `complex` type. It is an `intent(out)` argument.
286+
287+
### Syntax
288+
289+
`call ` [[stdlib_sparse_conversion(module):coo2dense(interface)]] `(coo,dense)`
290+
291+
### Arguments
292+
293+
`coo` : Shall be a `COO` type of `real` or `complex` type. It is an `intent(in)` argument.
294+
295+
`dense` : Shall be a rank-2 array of `real` or `complex` type. It is an `intent(out)` argument.
296+
297+
### Syntax
298+
299+
`call ` [[stdlib_sparse_conversion(module):coo2csr(interface)]] `(coo,csr)`
300+
301+
### Arguments
302+
303+
`coo` : Shall be a `COO` type of `real` or `complex` type. It is an `intent(in)` argument.
304+
305+
`csr` : Shall be a `CSR` type of `real` or `complex` type. It is an `intent(out)` argument.
306+
307+
### Syntax
308+
309+
`call ` [[stdlib_sparse_conversion(module):coo2csc(interface)]] `(coo,csc)`
310+
311+
### Arguments
312+
313+
`coo` : Shall be a `COO` type of `real` or `complex` type. It is an `intent(in)` argument.
314+
315+
`csc` : Shall be a `CSC` type of `real` or `complex` type. It is an `intent(out)` argument.
316+
317+
### Syntax
318+
319+
`call ` [[stdlib_sparse_conversion(module):csr2coo(interface)]] `(csr,coo)`
320+
321+
### Arguments
322+
323+
`csr` : Shall be a `CSR` type of `real` or `complex` type. It is an `intent(in)` argument.
324+
325+
`coo` : Shall be a `COO` type of `real` or `complex` type. It is an `intent(out)` argument.
326+
327+
### Syntax
328+
329+
`call ` [[stdlib_sparse_conversion(module):csr2sellc(interface)]] `(csr,sellc[,chunk])`
330+
331+
### Arguments
332+
333+
`csr` : Shall be a `CSR` type of `real` or `complex` type. It is an `intent(in)` argument.
334+
335+
`sellc` : Shall be a `SELLC` type of `real` or `complex` type. It is an `intent(out)` argument.
336+
337+
`chunk`, `optional`: chunk size for the Sliced ELLPACK format. It is an `intent(in)` argument.
338+
339+
### Syntax
340+
341+
`call ` [[stdlib_sparse_conversion(module):csr2sellc(interface)]] `(csr,ell[,num_nz_rows])`
342+
343+
### Arguments
344+
345+
`csr` : Shall be a `CSR` type of `real` or `complex` type. It is an `intent(in)` argument.
346+
347+
`ell` : Shall be a `ELL` type of `real` or `complex` type. It is an `intent(out)` argument.
348+
349+
`num_nz_rows`, `optional`: number of non zeros per row. If not give, it will correspond to the size of the longest row in the `CSR` matrix. It is an `intent(in)` argument.
350+
351+
### Syntax
352+
353+
`call ` [[stdlib_sparse_conversion(module):csc2coo(interface)]] `(csc,coo)`
354+
355+
### Arguments
356+
357+
`csc` : Shall be a `CSC` type of `real` or `complex` type. It is an `intent(in)` argument.
358+
359+
`coo` : Shall be a `COO` type of `real` or `complex` type. It is an `intent(out)` argument.
360+
361+
## Example
362+
```fortran
363+
{!example/linalg/example_sparse_spmv.f90!}
364+
```

‎example/linalg/CMakeLists.txt

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,9 @@ ADD_EXAMPLE(get_norm)
3333
ADD_EXAMPLE(solve1)
3434
ADD_EXAMPLE(solve2)
3535
ADD_EXAMPLE(solve3)
36+
ADD_EXAMPLE(sparse_from_ijv)
37+
ADD_EXAMPLE(sparse_data_accessors)
38+
ADD_EXAMPLE(sparse_spmv)
3639
ADD_EXAMPLE(svd)
3740
ADD_EXAMPLE(svdvals)
3841
ADD_EXAMPLE(determinant)
Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
program example_sparse_data_accessors
2+
use stdlib_linalg_constants, only: dp
3+
use stdlib_sparse
4+
implicit none
5+
6+
real(dp) :: mat(2,2)
7+
real(dp), allocatable :: dense(:,:)
8+
type(CSR_dp_type) :: CSR
9+
type(COO_dp_type) :: COO
10+
integer :: i, j, locdof(2)
11+
12+
! Initial data
13+
mat(:,1) = [1._dp,2._dp]
14+
mat(:,2) = [2._dp,1._dp]
15+
allocate(dense(5,5) , source = 0._dp)
16+
do i = 0, 3
17+
dense(1+i:2+i,1+i:2+i) = dense(1+i:2+i,1+i:2+i) + mat
18+
end do
19+
20+
print *, 'Original Matrix'
21+
do j = 1 , 5
22+
print '(5f8.1)',dense(j,:)
23+
end do
24+
25+
! Initialize CSR data and reset dense reference matrix
26+
call dense2coo(dense,COO)
27+
call coo2csr(COO,CSR)
28+
CSR%data = 0._dp
29+
dense = 0._dp
30+
31+
! Iteratively add blocks of data
32+
do i = 0, 3
33+
locdof(1:2) = [1+i,2+i]
34+
call CSR%add(locdof,locdof,mat)
35+
! lets print a dense view of every step
36+
call csr2dense(CSR,dense)
37+
print '(A,I2)', 'Add block :', i+1
38+
do j = 1 , 5
39+
print '(5f8.1)',dense(j,:)
40+
end do
41+
end do
42+
43+
! Request values from the matrix
44+
print *, ''
45+
print *, 'within sparse pattern :',CSR%at(2,1)
46+
print *, 'outside sparse pattern :',CSR%at(5,2)
47+
print *, 'outside matrix pattern :',CSR%at(7,7)
48+
49+
end program example_sparse_data_accessors
Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
program example_sparse_from_ijv
2+
use stdlib_linalg_constants, only: dp
3+
use stdlib_sparse
4+
implicit none
5+
6+
integer :: row(10), col(10)
7+
real(dp) :: data(10)
8+
type(COO_dp_type) :: COO
9+
type(CSR_dp_type) :: CSR
10+
type(ELL_dp_type) :: ELL
11+
integer :: i, j
12+
13+
! Initial data
14+
row = [1,1,2,2,3,3,3,4,4,4]
15+
col = [1,5,1,2,2,3,4,1,3,4]
16+
data = real([9,-3,4,7,8,-1,8,4,5,6] , kind = dp )
17+
18+
! Create a COO matrix from triplet
19+
call from_ijv(COO,row,col,data)
20+
print *, 'COO'
21+
print *, ' i, j, v'
22+
do i = 1, COO%nnz
23+
print '(2I4,f8.1)', COO%index(:,i), COO%data(i)
24+
end do
25+
26+
! Create a CSR matrix from triplet
27+
call from_ijv(CSR,row,col,data)
28+
print *, 'CSR'
29+
print '(A,5I8)', 'rowptr :', CSR%rowptr
30+
print '(A,10I8)', 'col :', CSR%col
31+
print '(A,10f8.1)', 'data :', CSR%data
32+
33+
! Create an ELL matrix from triplet
34+
call from_ijv(ELL,row,col,data)
35+
print *, 'ELL'
36+
print *, ' index | data'
37+
do i = 1, ELL%nrows
38+
print '(3I4,x,3f8.1)', ELL%index(i,:) , ELL%data(i,:)
39+
end do
40+
41+
end program example_sparse_from_ijv
Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
program example_sparse_spmv
2+
use stdlib_linalg_constants, only: dp
3+
use stdlib_sparse
4+
implicit none
5+
6+
integer, parameter :: m = 4, n = 2
7+
real(dp) :: A(m,n), x(n)
8+
real(dp) :: y_dense(m), y_coo(m), y_csr(m)
9+
real(dp) :: alpha, beta
10+
type(COO_dp_type) :: COO
11+
type(CSR_dp_type) :: CSR
12+
13+
call random_number(A)
14+
! Convert from dense to COO and CSR matrices
15+
call dense2coo( A , COO )
16+
call coo2csr( COO , CSR )
17+
18+
! Initialize vectors
19+
x = 1._dp
20+
y_dense = 2._dp
21+
y_coo = y_dense
22+
y_csr = y_dense
23+
24+
! Perform matrix-vector product
25+
alpha = 3._dp; beta = 2._dp
26+
y_dense = alpha * matmul(A,x) + beta * y_dense
27+
call spmv( COO , x , y_coo , alpha = alpha, beta = beta )
28+
call spmv( CSR , x , y_csr , alpha = alpha, beta = beta )
29+
30+
print *, 'dense :', y_dense
31+
print *, 'coo :', y_coo
32+
print *, 'csr :', y_csr
33+
34+
end program example_sparse_spmv

‎include/common.fypp

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@
3838

3939
#! Real types to be considered during templating
4040
#:set REAL_TYPES = ["real({})".format(k) for k in REAL_KINDS]
41+
#:set REAL_SUFFIX = REAL_KINDS
4142

4243
#! Collected (kind, type) tuples for real types
4344
#:set REAL_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_INIT))
@@ -75,6 +76,7 @@ $:"s" if cmplx=="c" else "d" if cmplx=="z" else "x" if cmplx=="y" else "q" if cm
7576

7677
#! Complex types to be considered during templating
7778
#:set CMPLX_TYPES = ["complex({})".format(k) for k in CMPLX_KINDS]
79+
#:set CMPLX_SUFFIX = ["c{}".format(k) for k in CMPLX_KINDS]
7880

7981
#! Collected (kind, type, initial) tuples for complex types
8082
#:set CMPLX_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_INIT))
@@ -115,6 +117,9 @@ $:"s" if cmplx=="c" else "d" if cmplx=="z" else "x" if cmplx=="y" else "q" if cm
115117
#! Bitset types to be considered during templating
116118
#:set BITSET_TYPES = ["type({})".format(k) for k in BITSET_KINDS]
117119

120+
#! Sparse types to be considered during templating
121+
#:set SPARSE_KINDS = ["COO", "CSR", "CSC", "ELL"]
122+
118123
#! Whether Fortran 90 compatible code should be generated
119124
#:set VERSION90 = defined('VERSION90')
120125

‎src/CMakeLists.txt

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -42,6 +42,10 @@ set(fppFiles
4242
stdlib_sorting_ord_sort.fypp
4343
stdlib_sorting_sort.fypp
4444
stdlib_sorting_sort_index.fypp
45+
stdlib_sparse_constants.fypp
46+
stdlib_sparse_conversion.fypp
47+
stdlib_sparse_kinds.fypp
48+
stdlib_sparse_spmv.fypp
4549
stdlib_specialfunctions_gamma.fypp
4650
stdlib_stats.fypp
4751
stdlib_stats_corr.fypp
@@ -114,6 +118,7 @@ set(SRC
114118
stdlib_logger.f90
115119
stdlib_sorting_radix_sort.f90
116120
stdlib_system.F90
121+
stdlib_sparse.f90
117122
stdlib_specialfunctions.f90
118123
stdlib_specialfunctions_legendre.f90
119124
stdlib_quadrature_gauss.f90

‎src/stdlib_sparse.f90

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
!! public API
2+
module stdlib_sparse
3+
use stdlib_sparse_kinds
4+
use stdlib_sparse_conversion
5+
use stdlib_sparse_spmv
6+
end module stdlib_sparse

‎src/stdlib_sparse_constants.fypp

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
#:include "common.fypp"
2+
#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
3+
#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
4+
module stdlib_sparse_constants
5+
use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp
6+
7+
implicit none
8+
public
9+
10+
enum, bind(C)
11+
enumerator :: sparse_full !! Full Sparse matrix (no symmetry considerations)
12+
enumerator :: sparse_lower !! Symmetric Sparse matrix with triangular inferior storage
13+
enumerator :: sparse_upper !! Symmetric Sparse matrix with triangular supperior storage
14+
end enum
15+
16+
character(1), parameter :: sparse_op_none = 'N' !! no transpose
17+
character(1), parameter :: sparse_op_transpose = 'T' !! transpose
18+
character(1), parameter :: sparse_op_hermitian = 'H' !! conjugate or hermitian transpose
19+
20+
! Integer size support for ILP64 builds should be done here
21+
integer, parameter :: ilp = int32
22+
23+
#:for k1, t1, s1 in (R_KINDS_TYPES)
24+
${t1}$, parameter :: zero_${s1}$ = 0._${k1}$
25+
${t1}$, parameter :: one_${s1}$ = 1._${k1}$
26+
#:endfor
27+
#:for k1, t1, s1 in (C_KINDS_TYPES)
28+
${t1}$, parameter :: zero_${s1}$ = (0._${k1}$,0._${k1}$)
29+
${t1}$, parameter :: one_${s1}$ = (1._${k1}$,1._${k1}$)
30+
#:endfor
31+
32+
end module stdlib_sparse_constants

‎src/stdlib_sparse_conversion.fypp

Lines changed: 938 additions & 0 deletions
Large diffs are not rendered by default.

‎src/stdlib_sparse_kinds.fypp

Lines changed: 593 additions & 0 deletions
Large diffs are not rendered by default.

‎src/stdlib_sparse_spmv.fypp

Lines changed: 593 additions & 0 deletions
Large diffs are not rendered by default.

‎test/linalg/CMakeLists.txt

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@ set(
1212
"test_linalg_qr.fypp"
1313
"test_linalg_svd.fypp"
1414
"test_linalg_matrix_property_checks.fypp"
15+
"test_linalg_sparse.fypp"
1516
)
1617
fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
1718

@@ -27,3 +28,4 @@ ADDTEST(linalg_lstsq)
2728
ADDTEST(linalg_qr)
2829
ADDTEST(linalg_svd)
2930
ADDTEST(blas_lapack)
31+
ADDTEST(linalg_sparse)

‎test/linalg/test_linalg_sparse.fypp

Lines changed: 414 additions & 0 deletions
Large diffs are not rendered by default.

0 commit comments

Comments
 (0)
Please sign in to comment.