Sparse Matrices¶
This chapter describes functions for the construction and manipulation of sparse matrices, matrices which are populated primarily with zeros and contain only a few non-zero elements. Sparse matrices often appear in the solution of partial differential equations. It is beneficial to use specialized data structures and algorithms for storing and working with sparse matrices, since dense matrix algorithms and structures can be very slow and use huge amounts of memory when applied to sparse matrices.
The header file gsl_spmatrix.h contains the prototypes for the
sparse matrix functions and related declarations.
Overview¶
These routines provide support for constructing and manipulating
sparse matrices in GSL, using an API similar to gsl_matrix.
The basic structure is called gsl_spmatrix. There are
three supported storage formats for sparse matrices: the triplet,
compressed column storage (CCS) and compressed row storage (CRS)
formats. The triplet format stores triplets
-
gsl_spmatrix¶ This structure is defined as:
typedef struct { size_t size1; size_t size2; size_t *i; double *data; size_t *p; size_t nzmax; size_t nz; gsl_spmatrix_tree *tree_data; void *work; size_t sptype; } gsl_spmatrix;
This defines a
, then and .size1-by-size2sparse matrix. The number of non-zero elements currently in the matrix is given bynz. For the triplet representation,i,p, anddataare arrays of sizenzwhich contain the row indices, column indices, and element value, respectively. So ifFor compressed column storage,
, then and .ianddataare arrays of sizenzcontaining the row indices and element values, identical to the triplet case.pis an array of sizesize2+ 1 wherep[j]points to the index indataof the start of columnj. Thus, ifFor compressed row storage,
, then and .ianddataare arrays of sizenzcontaining the column indices and element values, identical to the triplet case.pis an array of sizesize1+ 1 wherep[i]points to the index indataof the start of rowi. Thus, ifThe parameter
tree_datais a binary tree structure used in the triplet representation, specifically a balanced AVL tree. This speeds up element searches and duplicate detection during the matrix assembly process. The parameterworkis additional workspace needed for various operations like converting from triplet to compressed storage.sptypeindicates the type of storage format being used (triplet, CCS or CRS).The compressed storage format defined above makes it very simple to interface with sophisticated external linear solver libraries which accept compressed storage input. The user can simply pass the arrays
i,p, anddataas the inputs to external libraries.
Allocation¶
The functions for allocating memory for a sparse matrix follow the style of
malloc() and free(). They also perform their own error checking. If
there is insufficient memory available to allocate a matrix then the functions
call the GSL error handler with an error code of GSL_ENOMEM in addition
to returning a null pointer.
-
gsl_spmatrix *
gsl_spmatrix_alloc(const size_t n1, const size_t n2)¶ This function allocates a sparse matrix of size
. The functionn1-by-n2and initializes it to all zeros. If the size of the matrix is not known at allocation time, bothn1andn2may be set to 1, and they will automatically grow as elements are added to the matrix. This function sets the matrix to the triplet representation, which is the easiest for adding and accessing matrix elements. This function tries to make a reasonable guess for the number of non-zero elements (nzmax) which will be added to the matrix by assuming a sparse density ofgsl_spmatrix_alloc_nzmax()can be used if this number is known more accurately. The workspace is of size .
-
gsl_spmatrix *
gsl_spmatrix_alloc_nzmax(const size_t n1, const size_t n2, const size_t nzmax, const size_t sptype)¶ This function allocates a sparse matrix of size
n1-by-n2and initializes it to all zeros. If the size of the matrix is not known at allocation time, bothn1andn2may be set to 1, and they will automatically grow as elements are added to the matrix. The parameternzmaxspecifies the maximum number of non-zero elements which will be added to the matrix. It does not need to be precisely known in advance, since storage space will automatically grow usinggsl_spmatrix_realloc()ifnzmaxis not large enough. Accurate knowledge of this parameter reduces the number of reallocation calls required. The parametersptypespecifies the storage format of the sparse matrix. Possible values are-
GSL_SPMATRIX_TRIPLET¶ This flag specifies triplet storage.
-
GSL_SPMATRIX_CCS¶ This flag specifies compressed column storage.
-
GSL_SPMATRIX_CRS¶ This flag specifies compressed row storage.
The allocated
.gsl_spmatrixstructure is of size-
-
int
gsl_spmatrix_realloc(const size_t nzmax, gsl_spmatrix * m)¶ This function reallocates the storage space for
mto accomodatenzmaxnon-zero elements. It is typically called internally bygsl_spmatrix_set()if the user wants to add more elements to the sparse matrix than the previously specifiednzmax.
-
void
gsl_spmatrix_free(gsl_spmatrix * m)¶ This function frees the memory associated with the sparse matrix
m.
Accessing Matrix Elements¶
-
double
gsl_spmatrix_get(const gsl_spmatrix * m, const size_t i, const size_t j)¶ This function returns element (
i,j) of the matrixm. The matrix may be in triplet or compressed format.
-
int
gsl_spmatrix_set(gsl_spmatrix * m, const size_t i, const size_t j, const double x)¶ This function sets element (
i,j) of the matrixmto the valuex. The matrix must be in triplet representation.
-
double *
gsl_spmatrix_ptr(gsl_spmatrix * m, const size_t i, const size_t j)¶ This function returns a pointer to the (
i,j) element of the matrixm. If the (i,j) element is not explicitly stored in the matrix, a null pointer is returned.
Initializing Matrix Elements¶
Since the sparse matrix format only stores the non-zero elements, it is automatically
initialized to zero upon allocation. The function gsl_spmatrix_set_zero() may
be used to re-initialize a matrix to zero after elements have been added to it.
-
int
gsl_spmatrix_set_zero(gsl_spmatrix * m)¶ This function sets (or resets) all the elements of the matrix
mto zero.
Reading and Writing Matrices¶
-
int
gsl_spmatrix_fwrite(FILE * stream, const gsl_spmatrix * m)¶ This function writes the elements of the matrix
mto the streamstreamin binary format. The return value is 0 for success andGSL_EFAILEDif there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
-
int
gsl_spmatrix_fread(FILE * stream, gsl_spmatrix * m)¶ This function reads into the matrix
mfrom the open streamstreamin binary format. The matrixmmust be preallocated with the correct storage format, dimensions and have a sufficiently largenzmaxin order to read in all matrix elements, otherwiseGSL_EBADLENis returned. The return value is 0 for success andGSL_EFAILEDif there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
-
int
gsl_spmatrix_fprintf(FILE * stream, const gsl_spmatrix * m, const char * format)¶ This function writes the elements of the matrix
mline-by-line to the streamstreamusing the format specifierformat, which should be one of the%g,%eor%fformats for floating point numbers. The function returns 0 for success andGSL_EFAILEDif there was a problem writing to the file. The input matrixmmay be in any storage format, and the output file will be written in MatrixMarket format.
-
gsl_spmatrix *
gsl_spmatrix_fscanf(FILE * stream)¶ This function reads sparse matrix data in the MatrixMarket format from the stream
streamand stores it in a newly allocated matrix which is returned in triplet format. The function returns 0 for success andGSL_EFAILEDif there was a problem reading from the file. The user should free the returned matrix when it is no longer needed.
Copying Matrices¶
-
int
gsl_spmatrix_memcpy(gsl_spmatrix * dest, const gsl_spmatrix * src)¶ This function copies the elements of the sparse matrix
srcintodest. The two matrices must have the same dimensions and be in the same storage format.
Exchanging Rows and Columns¶
-
int
gsl_spmatrix_transpose_memcpy(gsl_spmatrix * dest, const gsl_spmatrix * src)¶ This function copies the transpose of the sparse matrix
srcintodest. The dimensions ofdestmust match the transpose of the matrixsrc. Also, both matrices must use the same sparse storage format.
-
int
gsl_spmatrix_transpose(gsl_spmatrix * m)¶ This function replaces the matrix
mby its transpose, preserving the storage format of the input matrix. Currently, only triplet matrix inputs are supported.
-
int
gsl_spmatrix_transpose2(gsl_spmatrix * m)¶ This function replaces the matrix
mby its transpose, but changes the storage format for compressed matrix inputs. Since compressed column storage is the transpose of compressed row storage, this function simply converts a CCS matrix to CRS and vice versa. This is the most efficient way to transpose a compressed storage matrix, but the user should note that the storage format of their compressed matrix will change on output. For triplet matrices, the output matrix is also in triplet storage.
Matrix Operations¶
-
int
gsl_spmatrix_add(gsl_spmatrix * c, const gsl_spmatrix * a, const gsl_spmatrix * b)¶ This function computes the sum
. The three matrices must have the same dimensions and be stored in a compressed format.
-
int
gsl_spmatrix_scale(gsl_spmatrix * m, const double x)¶ This function scales all elements of the matrix
is stored inmby the constant factorx. The resultm.
Matrix Properties¶
-
size_t
gsl_spmatrix_nnz(const gsl_spmatrix * m)¶ This function returns the number of non-zero elements in
m.
-
int
gsl_spmatrix_equal(const gsl_spmatrix * a, const gsl_spmatrix * b)¶ This function returns 1 if the matrices
aandbare equal (by comparison of element values) and 0 otherwise. The matricesaandbmust be in the same sparse storage format for comparison.
Finding Maximum and Minimum Elements¶
-
int
gsl_spmatrix_minmax(const gsl_spmatrix * m, double * min_out, double * max_out)¶ This function returns the minimum and maximum elements of the matrix
m, storing them inmin_outandmax_out, and searching only the non-zero values.
Compressed Format¶
GSL supports compressed column storage (CCS) and compressed row storage (CRS) formats.
-
gsl_spmatrix *
gsl_spmatrix_ccs(const gsl_spmatrix * T)¶ This function creates a sparse matrix in compressed column format from the input sparse matrix
Twhich must be in triplet format. A pointer to a newly allocated matrix is returned. The calling function should free the newly allocated matrix when it is no longer needed.
-
gsl_spmatrix *
gsl_spmatrix_crs(const gsl_spmatrix * T)¶ This function creates a sparse matrix in compressed row format from the input sparse matrix
Twhich must be in triplet format. A pointer to a newly allocated matrix is returned. The calling function should free the newly allocated matrix when it is no longer needed.
Conversion Between Sparse and Dense Matrices¶
The gsl_spmatrix structure can be converted into the dense gsl_matrix
format and vice versa with the following routines.
-
int
gsl_spmatrix_d2sp(gsl_spmatrix * S, const gsl_matrix * A)¶ This function converts the dense matrix
Ainto sparse triplet format and stores the result inS.
-
int
gsl_spmatrix_sp2d(gsl_matrix * A, const gsl_spmatrix * S)¶ This function converts the sparse matrix
Sinto a dense matrix and stores the result inA.Smust be in triplet format.
Examples¶
The following example program builds a 5-by-4 sparse matrix and prints it in triplet, compressed column, and compressed row format. The matrix which is constructed is
The output of the program is:
printing all matrix elements:
A(0,0) = 0
A(0,1) = 0
A(0,2) = 3.1
A(0,3) = 4.6
A(1,0) = 1
.
.
.
A(4,0) = 4.1
A(4,1) = 0
A(4,2) = 0
A(4,3) = 0
matrix in triplet format (i,j,Aij):
(0, 2, 3.1)
(0, 3, 4.6)
(1, 0, 1.0)
(1, 2, 7.2)
(3, 0, 2.1)
(3, 1, 2.9)
(3, 3, 8.5)
(4, 0, 4.1)
matrix in compressed column format:
i = [ 1, 3, 4, 3, 0, 1, 0, 3, ]
p = [ 0, 3, 4, 6, 8, ]
d = [ 1, 2.1, 4.1, 2.9, 3.1, 7.2, 4.6, 8.5, ]
matrix in compressed row format:
i = [ 2, 3, 0, 2, 0, 1, 3, 0, ]
p = [ 0, 2, 4, 4, 7, 8, ]
d = [ 3.1, 4.6, 1, 7.2, 2.1, 2.9, 8.5, 4.1, ]
We see in the compressed column output, the data array stores each column contiguously, the array
stores the row index of the corresponding data element, and the array stores the index of the start of each column in the data array. Similarly, for the compressed row output, the data array stores each row contiguously, the array stores the column index of the corresponding data element, and the array stores the index of the start of each row in the data array.#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_spmatrix.h>
int
main()
{
gsl_spmatrix *A = gsl_spmatrix_alloc(5, 4); /* triplet format */
gsl_spmatrix *B, *C;
size_t i, j;
/* build the sparse matrix */
gsl_spmatrix_set(A, 0, 2, 3.1);
gsl_spmatrix_set(A, 0, 3, 4.6);
gsl_spmatrix_set(A, 1, 0, 1.0);
gsl_spmatrix_set(A, 1, 2, 7.2);
gsl_spmatrix_set(A, 3, 0, 2.1);
gsl_spmatrix_set(A, 3, 1, 2.9);
gsl_spmatrix_set(A, 3, 3, 8.5);
gsl_spmatrix_set(A, 4, 0, 4.1);
printf("printing all matrix elements:\n");
for (i = 0; i < 5; ++i)
for (j = 0; j < 4; ++j)
printf("A(%zu,%zu) = %g\n", i, j,
gsl_spmatrix_get(A, i, j));
/* print out elements in triplet format */
printf("matrix in triplet format (i,j,Aij):\n");
gsl_spmatrix_fprintf(stdout, A, "%.1f");
/* convert to compressed column format */
B = gsl_spmatrix_ccs(A);
printf("matrix in compressed column format:\n");
printf("i = [ ");
for (i = 0; i < B->nz; ++i)
printf("%zu, ", B->i[i]);
printf("]\n");
printf("p = [ ");
for (i = 0; i < B->size2 + 1; ++i)
printf("%zu, ", B->p[i]);
printf("]\n");
printf("d = [ ");
for (i = 0; i < B->nz; ++i)
printf("%g, ", B->data[i]);
printf("]\n");
/* convert to compressed row format */
C = gsl_spmatrix_crs(A);
printf("matrix in compressed row format:\n");
printf("i = [ ");
for (i = 0; i < C->nz; ++i)
printf("%zu, ", C->i[i]);
printf("]\n");
printf("p = [ ");
for (i = 0; i < C->size1 + 1; ++i)
printf("%zu, ", C->p[i]);
printf("]\n");
printf("d = [ ");
for (i = 0; i < C->nz; ++i)
printf("%g, ", C->data[i]);
printf("]\n");
gsl_spmatrix_free(A);
gsl_spmatrix_free(B);
gsl_spmatrix_free(C);
return 0;
}
References and Further Reading¶
The algorithms used by these functions are described in the following sources,
- Davis, T. A., Direct Methods for Sparse Linear Systems, SIAM, 2006.
- CSparse software library, https://www.cise.ufl.edu/research/sparse/CSparse