Sparse matrix¶
This library implements a lightweight class smatrix
for the manipulation of sparse matrices, using solely the castor::matrix
class without other external libraries. Many builders, functions and operators are available in a similar fashion as for the Class matrix, see the smatrix
API. Most of the manipulations should be transparent to the user and the feeling should be really Matlab-like. Basic examples are given in the corresponding Examples section.
The smatrix
class is installed automatically with the other headers of castor library.
This is work-in-progress
Examples¶
In order to use a smatrix
, the corresponding header castor/smatrix.hpp
needs to be included. A minimum working example would thus look like this:
#include "castor/matrix.hpp"
#include "castor/smatrix.hpp"
int main()
{
smatrix<> s;
// your code here
return 0;
}
Building a smatrix¶
A smatrix
will mostly behave like a normal dense matrix
. For example, let us create a 4 x 4
smatrix
of double
.
smatrix<> As;
disp(As,2);
Sparse matrix 4x4 of type 'd' with 0 elements (0 B):
-empty-
The newly created smatrix
is empty. It is rather easy to set the values of the entries. For example:
As(1,2) = -0.5;
disp(As,2);
Sparse matrix 4x4 of type 'd' with 1 elements (8 B):
(1,2) -0.5
Now As
contains one non-zero element with the bilinear index (1,2)
(second line, third column). However, this way of filling a smatrix
is definitely not recommended as it involves a lot of memory management and may affect dramatically the performances. One way to do so is to first create the matrix in triplet format (I,J,VALUE)
and only afterward build the corresponding smatrix
as illustrated below:
matrix<std::size_t> I = {0,2,3};
matrix<std::size_t> J = {1,1,2};
matrix<> V = {0.5, -1/3., M_PI};
As = smatrix<>(I,J,V,4,4);
disp(As,2);
Sparse matrix 4x4 of type 'd' with 3 elements (24 B):
(0,1) 0.5
(2,1) -0.333333
(3,2) 3.14159
Yes, we reaffected As
to a new smatrix
. The old data is automatically discarded so one should be careful when performing such an action. As for matrix
, it is possible to clear the content of a smatrix
(the object is reinitialized). Let us now add a 0
.
As(3,3) = 0.; // but, why ?
disp(As,2);
Sparse matrix 4x4 of type 'd' with 4 elements (32 B):
(0,1) 0.5
(2,1) -0.333333
(3,2) 3.14159
(3,3) 0
A zero value is added. In order to clean a smatrix
, a simple call to check
is sufficient:
check(As);
disp(As,2);
Sparse matrix 4x4 of type 'd' with 3 elements (24 B):
(0,1) 0.5
(2,1) -0.333333
(3,2) 3.14159
Everything went back to normal! Now, let us use one of the builders in order to create an identity sparse matrix. It is also possible to convert back to the triplet format.
auto Bs = speye<>(4,4);
disp(Bs,2);
matrix<std::size_t> IB,JB;
matrix<> VB;
std::tie(IB,JB,VB) = find(Bs);
disp(transpose(vertcat(vertcat(IB,JB),VB)),2);
Sparse matrix 4x4 of type 'd' with 4 elements (32 B):
(0,0) 1
(1,1) 1
(2,2) 1
(3,3) 1
Matrix 4x3 of type 'd' (96 B):
0 0 1.00000
1.00000 1.00000 1.00000
2.00000 2.00000 1.00000
3.00000 3.00000 1.00000
The matrices IB,JB,VB
are returned as line vectors. To obtain a better display, we concatenated them vertically and tranposed the result.
Basic manipulations¶
In this section, we start with start from scratch so everything written in the previous section should be discarded from your main
function. Let us create two matrices
smatrix<> As = speye(4,4);
matrix<std::size_t> I({1,1,2,2}), J({1,2,1,2});
matrix<> V({1.,1.,1.,1.});
smatrix<> Bs = smatrix<>(I,J,V,4,4);
As
is a 4 x 4
identity matrix and Bs
is a matrix with the interior filled with ones. Here is an example of basic manipulations:
auto Cs = 1.5*As - Bs/2.;
disp(Cs,2);
Sparse matrix 4x4 of type 'd' with 6 elements (48 B):
(0,0) 1.5
(1,1) 1
(1,2) -0.5
(2,1) -0.5
(2,2) 1
(3,3) 1.5
What is the number of non-zero elements, again ?
std::cout << "nnz(Cs) = " << nnz(Cs) << std::endl;
nnz(Cs) = 6
It is possible to get the value of any entry:
std::cout << "Cs(1,2) = " << Cs(1,2) << std::endl;
std::cout << "Cs(1,3) = " << Cs(1,3) << std::endl;
Cs(1,2) = -0.5
Cs(1,3) = 0
Now, let us multiply Cs
by a 4 x 4
dense matrix
:
auto D = mtimes(Cs,ones<>(4));
disp(D,2); // :)
Matrix 4x4 of type 'd' (128 B):
1.50000 1.50000 1.50000 1.50000
0.50000 0.50000 0.50000 0.50000
0.50000 0.50000 0.50000 0.50000
1.50000 1.50000 1.50000 1.50000
One last manipulation and we are good for this example.
Cs(0,3) = M_PI;
auto Es = Cs - transpose(Cs);
check(Es); // drop the zeros
disp(Es,2);
disp(full(Es,2));
Sparse matrix 4x4 of type 'd' with 2 elements (16 B):
(0,3) 3.14159
(3,0) -3.14159
Matrix 4x4 of type 'd' (128 B):
0 0 0 3.14159
0 0 0 0
0 0 0 0
-3.14159 0 0 0