You are on page 1of 28

# Benchmark of C++ Libraries for Sparse Matrix

Computation
Georg Holzmann
http://grh.mur.at
email: grh@mur.at
August 2007

This report presents benchmarks of C++ scientific computing libraries for small
and medium size sparse matrices. Be warned: these benchmarks are very special-
ized on a neural network like algorithm I had to implement. However, also the
initialization time of sparse matrices and a matrix-vector multiplication was mea-
sured, which might be of general interest.

WARNING
At the time of its writing this document did not consider the eigen
library (http://eigen.tuxfamily.org).
Please evaluate also this very nice, fast and well maintained library
before making any decisions!
information.

Contents
1 Introduction 2

2 Benchmarks 2
2.1 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Matrix-Vector Multiplication . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 Neural Network like Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3 Results 4

4 Conclusion 12

## A Libraries and Flags 12

B Implementations 13

1
1 Introduction
Quite a lot open-source libraries exist for scientific computing, which makes it hard to choose
between them. Therefore, after communication on various mailing lists, I decided to perform
some benchmarks, specialized to the kind of algorithms I had to implement.
Ideally algorithms should be implemented in an object oriented, reuseable way, but of course
without a loss of performance. BLAS (Basic Linear Algebra Subprograms - see ) and LA-
PACK (Linear Algebra Package - see ) are the standard building blocks for efficient scientific
software. Highly optimized BLAS implementations for all relevant hardware architectures are
available, traditionally expressed in Fortran routines for scalar, vector and matrix operations
(called Level 1, 2 and 3 BLAS).
However, some of the here presented libraries still interface to BLAS without loosing the ob-
ject oriented concepts of component reusability, general programming, etc. Good performance
can be achieved with expression templates, closures and other design patterns (see  and ).
I have to note that these benchmarks are very special to the algorithms I had to implement
(neural network like operations) and should therefore not be seen as a general benchmark of
these libraries. In particular I benchmarked a (sparse) matrix-vector multiplication (see 2.2),
a neural network like operation (see 2.3) and the initialization time of the sparse matrix (see
2.1). Results to these examples can be found in section 3. A more general benchmark of
non-sparse algorithms can be found at http://projects.opencascade.org/btl/.

2 Benchmarks
I used the C(++) function clock() from ctime.h for the benchmarks (see ). As pointed out
in  the disadvantage with clock() is that the resolution is pretty low. Therefore one has to
repeat the code many times, where the data can stay in the cache and so the code runs faster.
However, I decided to use that method because in my application I also have to repeat the
operations in a similar loop and therefore the resolution was high enough.
The benchmarks were made with floats and doubles with sparse and nonsparse versions of
each library, where the term “connectivity” always specifies how many percent of the matrix
are non-zero elements.
The y-axis in the plots is usually given in CPU ticks per element. This means the number
of clock ticks elapsed between two clock() calls divided through the matrix size N (e.g. divided
through 100 for a 100x100 matrix). If one would want to calculate the absolute time needed
by an operation (in sec) one could use the following formula:

CP U T icksP erElement ∗ N
AbsT ime =
CLOCKS P ER SEC
where CLOCKS PER SEC is a macro constant expression of ctime.h, which specifies the re-
lation between a clock tick and a second (on the benchmark system: CLOCKS PER SEC=1000000).
As already said, each operation was repeated in a loop for one benchmark (see Table 1 for
the exact values). Additionally the whole benchmarks were repeated six times and averaged.
Table 2 shows some information about the benchmark system.
All libraries were compiled without debug mode with g++ (GCC) 4.1.2. The detailed
compiler and linker flags can be found in Appendix A, the source code in Appendix B.

2
Matrix Size (N) Nr. of repetitions
10 200000
30 100000
50 100000
70 50000
85 50000
100 10000
200 10000
300 6000
500 3000
1000 1000
2000 250
5000 50

## Table 1: Number of repetitions for a matrix size M.

OS linux 2.6.20
RAM 512 MB
Processor Intel(R) Pentium(R) 4 CPU 3.06GHz
CPU MHz 3056.931
Cache Size 512 KB

## Table 2: Benchmark system information.

2.1 Initialization
In the initialization phase random values are assigned to a matrix and three vectors with a
specific connectivity, like in the following code example:
for ( int i =0; i <N; ++i )
{
// v e c t o r s
i n ( i ) = ( r a n d v a l ()< i c o n n ) ? r a n d v a l ( ) : 0 ;
back ( i ) = ( r a n d v a l ()< f b c o n n ) ? r a n d v a l ( ) : 0 ;
x( i ) = randval ( ) ;

// m a t r i x
for ( int j =0; j <N; ++j )
{
W( i , j ) = ( r a n d v a l ()< conn ) ? r a n d v a l ( ) : 0 ;
}
}
This is mainly interesting for sparse matrix libraries, because they can need quite a lot of
time in the initialization phase (which might be critical or not).

3
2.2 Matrix-Vector Multiplication
The following matrix-vector product was benchmarked:

x=W ∗x
where W is a NxN sparse matrix with a specific connectivity and x is a Nx1 dense vector.
A temporary object is needed for such an operation, which results in an additional copy of the
vector x:

t=x
x=W ∗t
In some libraries the temporary object is generated automatically (which is the usual C++
way, but can result in a big loss of performance), some give warnings and some even produce
compilation errors and one has to code the temporary object explicitly.

## 2.3 Neural Network like Operation

The main performance critical loop of the neural network I had to implement (an Echo State
Network) looks like this:

x = α ∗ in + W ∗ x + β ∗ back
where in and back are vectors of size Nx1 and α, β are scalars. Now a kernel function is
applied to all the elements of vector x:

x = tanh(x)
And finally a dot product is calculated between vector out of size 2Nx1 and two concatenated
x vectors, where the last one is squared:

y = dot(out, (x, x2 ))

3 Results
In this section I want to show some illustrative results. For plots with all the tested parameters
and higher resolution see http://grh.mur.at/misc/sparse_benchmark_plots/.
Initialization benchmarks can be found in Figure 1 and 2. In general sparse matrices need of
course much more time than their nonsparse equivalents: in my examples FLENS took most
CPU time for initialization, followed by MTL and then Gmm++.
The matrix-vector product, with doubles and floats, can be seen in Figure 3, 4, 5, 6, 7 and
8. Usually Gmm++ is a little bit faster than FLENS, followed by MTL. Newmat and Blitz++
performed worse than the rest.
Results for the neural network like operation are shown in Figure 9, 10, 11, 12, 13 and
14. They are quite similar to the matrix-vector product, because this is the most expensive
operation. However, for small connectivities FLENS outperforms Gmm++, MTL is also quite
similar to Gmm++.

4
Figure 1: Initialization with doubles, connectivity 10%.

## Figure 2: Initialization with doubles, connectivity 40%.

5
Figure 3: matrix-vector product with doubles, connectivity 10%. NOTE: logarithmic y axis !

Figure 4: matrix-vector product with doubles, connectivity 10%, zoomed to lower area.

6
Figure 5: matrix-vector product with floats, connectivity 10%, zoomed to lower area.

Figure 6: matrix-vector product with doubles, connectivity 40%. NOTE: logarithmic y axis !

7
Figure 7: matrix-vector product with doubles, connectivity 40%, zoomed to lower area.

Figure 8: matrix-vector product with floats, connectivity 40%, zoomed to lower area.

8
Figure 9: NN like operation with doubles, connectivity 10%. NOTE: logarithmic y axis !

Figure 10: NN like operation with doubles, connectivity 10%, zoomed to lower area.

9
Figure 11: NN like operation with floats, connectivity 10%, zoomed to lower area.

Figure 12: NN like operation with doubles, connectivity 40%. NOTE: logarithmic y axis !

10
Figure 13: NN like operation with doubles, connectivity 40%, zoomed to lower area.

Figure 14: NN like operation with floats, connectivity 40%, zoomed to lower area.

11
4 Conclusion
In general one can say that Gmm++, FLENS and MTL have all very good performance.
Also the documentation is quite good for all these libraries. However, FLENS can be used in a
more intuitive way (not always), because it overloads operators +, *, -, ... without performance
overhead. Gmm++ was in my humble opinion easier to use than MTL, although the commands
are quite similar, but Gmm++ supports also sub-matrices for any matrix type.
Both Gmm++ and FLENS have an interface to BLAS and LAPACK, so it’s possible to use
vendor tuned libraries. In FLENS the bindings are activated per default, in Gmm++ one has
to define them. I tried linking Gmm++ against BLAS and LAPACK, but got no performance
boost also with non-sparse operations. FLENS was linked against ATLAS BLAS 1 . I tried also
GotoBLAS 2 , which is said to be faster than ATLAS, but as these BLAS implementations
usually come with no sparse matrix algorithms they perform worse.
According to the benchmarks for non-sparse algorithms presented at http://projects.
opencascade.org/btl/ MTL is even faster than the Fortran BLAS implementation and than
vendor tuned libraries like the Intel MKL 3 . Intel MKL supports also sparse matrices, which I
tried to test, but I did not get a free copy of it (because of some strange errors).

## A Libraries and Flags

The following libraries were tested:

## • Blitz++, see http://www.oonumerics.org/blitz/

• boost uBLAS, see www.boost.org/libs/numeric/
• FLENS, see http://flens.sourceforge.net/
• Gmm++, see http://home.gna.org/getfem/gmm_intro
• Matrix Template Library (MTL), see http://www.osl.iu.edu/research/mtl/
• newmat, see http://www.robertnz.net/nm_intro.htm

## Compilation/Linking Flags were set to:

• Blitz++:
CPPFLAGS = -msse -march=pentium4 -mfpmath=sse -funroll-loops -O3 -ftemplate-
depth-50;
LDFLAGS = -lc -lblitz -lm;
• boost uBLAS: CPPFLAGS = -DNDEBUG -ftemplate-depth-50 -msse -march=pentium4
-mfpmath=sse -funroll-loops -O3;
LDFLAGS 4 = -lc;
1
ATLAS: Automatic Tuned Linear Algebra Software, a portable, efficient BLAS implementation; http://
math-atlas.sourceforge.net/
2
GotoBLAS: highly optimized BLAS implementation; http://www.tacc.utexas.edu/resources/software/
#blas
3
Intel Math Kernel Library: http://www.intel.com/cd/software/products/asmo-na/eng/307757.htm
4
It did nont try linking boost uBLAS to BLAS and LAPACK. There exist such bindings but I guess especially
for sparse matrices they do not perform better, it should be tried though.

12
• FLENS:
CPPFLAGS = -g -O3 -fPIC -DNDEBUG -msse -march=pentium4 -mfpmath=sse -
funroll-loops -ftemplate-depth-50;
LDFLAGS = -llapack -latlas -lblas;

• Gmm++:
CPPFLAGS = -DNDEBUG -ftemplate-depth-50 -msse -march=pentium4 -mfpmath=sse
-funroll-loops -O3
LDFLAGS 5 = -lc

## • MTL: CPPFLAGS = -DNDEBUG -ftemplate-depth-50 -msse -march=pentium4 -mfpmath=sse

-funroll-loops -O3
LDFLAGS = -lc

## • newmat: CPPFLAGS = -ftemplate-depth-50 -msse -march=pentium4 -mfpmath=sse -

funroll-loops -O3
LDFLAGS = -lnewmat

B Implementations
Here the exact implementations of the test algorithms in all libraries are shown. Don’t be
confused about the little bit cryptic pointer notation!
benchmarks.tar.gz.

## Listing 1: Blitz++ Non-Sparse Implementation

// BLITZ++ NONSPARSE IMPLEMENTATION

## #include < b l i t z / a r r a y . h>

#include < b l i t z / t i n y v e c . h>

#include ” . . / bench−t o o l s . h”

using namespace b l i t z ;

## typedef Array<TYPE,2 > DMatrix ;

typedef Array<TYPE,1 > DVector ;

DMatrix ∗W;
DVector ∗ in ;
DVector ∗ back ;
DVector ∗ out ;
DVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
5
Here Gmm++ was not linked against BLAS and LAPACK. I tried this but the performance was worse or the
same for small/medium matrices.

13
{
W = new DMatrix (N,N ) ;
i n = new DVector (N ) ;
back = new DVector (N ) ;
out = new DVector ( 2 ∗N ) ;
x = new DVector (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
( ∗ i n ) ( i ) = ( r a n d v a l () < i c o n n ) ? r a n d v a l ( ) : 0 ;
( ∗ back ) ( i ) = ( r a n d v a l () < f b c o n n ) ? r a n d v a l ( ) : 0 ;
(∗ x ) ( i ) = randval ( ) ;

## f o r ( i n t j =0; j <N; ++j )

{
( ∗W) ( i , j ) = ( r a n d v a l () < conn ) ? r a n d v a l ( ) : 0 ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) ( i ) = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )
{
firstIndex i ;
secondIndex j ;

while ( i t e r −−)
{
// x = W ∗ x i s q u i t e t r i c k y
∗x = sum ( ( ∗W) ( j , i ) ∗ ( ∗ x ) ( i ) , j ) ;
// t h e r e s t i s e a s y
∗x += a l p h a ∗ ∗ i n ;
∗x += b e t a ∗ ∗ back ;

## f o r ( i n t i =0; i <N; ++i )

( ∗ x ) ( i ) = tanh ( ( ∗ x ) ( i ) ) ;

## // make a temporary v e c t o r ( x , x ˆ2)

y = sum ( ( ∗ out ) ( Range ( 0 ,N−1)) ∗ ∗x ) ;
f o r ( i n t i =0; i <N; ++i )
y += ( ∗ out ) (N+i ) ∗ pow ( ( ∗ x ) ( i ) , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
firstIndex i ;
secondIndex j ;

while ( i t e r −−)
{
// x = W ∗ x was much s l o w e r as a l l t h e compute ( )
// f u n c t i o n −> why ?

14
// so now I c a l c i n = W ∗ x

∗ i n = sum ( ( ∗W) ( j , i ) ∗ ( ∗ x ) ( i ) , j ) ;
}
}

## Listing 2: Boost Non-Sparse Implementation

// BOOST NONSPARSE IMPLEMENTATION

## #include <b o o s t / numeric / u b l a s / v e c t o r . hpp>

#include <b o o s t / numeric / u b l a s / matrix . hpp>
#include <b o o s t / numeric / u b l a s / i o . hpp>
#include <b o o s t / numeric / u b l a s / v e c t o r p r o x y . hpp>

#include ” . . / bench−t o o l s . h”

## typedef matrix<TYPE> GEMatrix ;

typedef v e c t o r <TYPE> DEVector ;

GEMatrix ∗W;
DEVector ∗ in ;
DEVector ∗ back ;
DEVector ∗ out ;
DEVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new GEMatrix (N,N ) ;
i n = new DEVector (N ) ;
back = new DEVector (N ) ;
out = new DEVector ( 2 ∗N ) ;
x = new DEVector (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
( ∗ i n ) ( i ) = ( r a n d v a l () < i c o n n ) ? r a n d v a l ( ) : 0 ;
( ∗ back ) ( i ) = ( r a n d v a l () < f b c o n n ) ? r a n d v a l ( ) : 0 ;
(∗ x ) ( i ) = randval ( ) ;

## f o r ( i n t j =0; j <N; ++j )

{
( ∗W) ( i , j ) = ( r a n d v a l () < conn ) ? r a n d v a l ( ) : 0 ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) ( i ) = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )

15
{
while ( i t e r −−)
{
∗x = prod ( ∗W, ∗ x ) ;
∗x += a l p h a ∗ ∗ i n ;
∗x += b e t a ∗ ∗ back ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =0; i <N; ++i )
( ∗ x ) ( i ) = tanh ( ( ∗ x ) ( i ) ) ;

// c a l c u l a t e y = [ o u t ] ∗ [ x ; x ˆ 2 ]
y = i n n e r p r o d ( s u b r a n g e ( ∗ out , 0 , N) , s u b r a n g e ( ∗ x , 0 , N) ) ;
f o r ( i n t i =0; i <N; ++i )
y += ( ∗ out ) (N+i ) ∗ pow ( ( ∗ x ) ( i ) , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
while ( i t e r −−)
{
∗x = prod ( ∗W, ∗ x ) ;
}
}

## Listing 3: Boost Sparse Implementation

// BOOST SPARSE IMPLEMENTATION

## #include <b o o s t / numeric / u b l a s / v e c t o r . hpp>

#include <b o o s t / numeric / u b l a s / v e c t o r s p a r s e . hpp>
#include <b o o s t / numeric / u b l a s / m a t r i x s p a r s e . hpp>
#include <b o o s t / numeric / u b l a s / i o . hpp>
#include <b o o s t / numeric / u b l a s / v e c t o r p r o x y . hpp>

#include ” . . / bench−t o o l s . h”

## typedef c o m p r e s s e d m a t r i x <TYPE> SMatrix ;

typedef c o m p r e s s e d v e c t o r <TYPE> SVector ;
typedef v e c t o r <TYPE> DEVector ;

SMatrix ∗W;
SVector ∗ i n ;
SVector ∗ back ;
DEVector ∗ out ;
DEVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )

16
{
W = new SMatrix (N,N ) ;
i n = new SVector (N ) ;
back = new SVector (N ) ;
out = new DEVector ( 2 ∗N ) ;
x = new DEVector (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
i f ( r a n d v a l () < i c o n n ) ( ∗ i n ) ( i ) = r a n d v a l ( ) ;
i f ( r a n d v a l () < f b c o n n ) ( ∗ back ) ( i ) = r a n d v a l ( ) ;
(∗ x ) ( i ) = randval ( ) ;

## f o r ( i n t j =0; j <N; ++j )

{
i f ( r a n d v a l () < conn ) ( ∗W) ( i , j ) = r a n d v a l ( ) ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) ( i ) = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )
{
while ( i t e r −−)
{
∗x = prod ( ∗W, ∗ x ) ;
∗x += a l p h a ∗ ∗ i n ;
∗x += b e t a ∗ ∗ back ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =0; i <N; ++i )
( ∗ x ) ( i ) = tanh ( ( ∗ x ) ( i ) ) ;

// c a l c u l a t e y = [ o u t ] ∗ [ x ; x ˆ 2 ]
y = i n n e r p r o d ( s u b r a n g e ( ∗ out , 0 , N) , s u b r a n g e ( ∗ x , 0 , N) ) ;
f o r ( i n t i =0; i <N; ++i )
y += ( ∗ out ) (N+i ) ∗ pow ( ( ∗ x ) ( i ) , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
while ( i t e r −−)
{
∗x = prod ( ∗W, ∗ x ) ;
}
}

## Listing 4: FLENS Non-Sparse Implementation

// FLENS NONSPARSE IMPLEMENTATION

## #include < f l e n s / f l e n s . h>

17
#include <math . h>
#include ” . . / bench−t o o l s . h”

using namespace f l e n s ;
using namespace s t d ;

## typedef GeMatrix<F u l l S t o r a g e <TYPE, ColMajor> > GEMatrix ;

typedef DenseVector<Array<TYPE> > DEVector ;

GEMatrix ∗W;
DEVector ∗ in ;
DEVector ∗ back ;
DEVector ∗ out ;
DEVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new GEMatrix (N,N ) ;
i n = new DEVector (N ) ;
back = new DEVector (N ) ;
out = new DEVector ( 2 ∗N ) ;
x = new DEVector (N ) ;

## f o r ( i n t i =1; i<=N; ++i )

{
i f ( r a n d v a l () < i c o n n ) ( ∗ i n ) ( i ) = r a n d v a l ( ) ;
i f ( r a n d v a l () < f b c o n n ) ( ∗ back ) ( i ) = r a n d v a l ( ) ;
(∗ x ) ( i ) = randval ( ) ;

## f o r ( i n t j =1; j<=N; ++j )

{
i f ( r a n d v a l () < conn ) ( ∗W) ( i , j ) = r a n d v a l ( ) ;
}
}

## f o r ( i n t i =1; i <=2∗N; ++i )

( ∗ out ) ( i ) = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )
{
// temporary o b j e c t s
DEVector t ;
DEVector : : View out1 = ( ∗ out ) ( ( 1 ,N) ) , out2 = ( ∗ out ) ( (N+1 ,2∗N ) ) ;

while ( i t e r −−)
{
t = ∗x ; // temporary o b j e c t needed f o r BLAS
∗x = a l p h a ∗ ∗ i n + ∗W ∗ t + b e t a ∗ ∗ back ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =1; i<=N; ++i )

18
( ∗ x ) ( i ) = tanh ( ( ∗ x ) ( i ) ) ;

// c a l c u l a t e y = [ o u t ] ∗ [ x ; x ˆ 2 ]
y = dot ( out1 , ∗x ) ;
f o r ( i n t i =1; i<=N; ++i )
y += out2 ( i ) ∗ pow ( ( ∗ x ) ( i ) , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
DEVector : : View t 1 = ( ∗ x ) ( ( 1 ,N ) ) ;

while ( i t e r −−)
{
∗x = ∗W ∗ t 1 ;
}
}

## Listing 5: FLENS Sparse Implementation

// FLENS SPARSE IMPLEMENTATION

## #include < f l e n s / f l e n s . h>

#include <math . h>
#include ” . . / bench−t o o l s . h”

using namespace f l e n s ;
using namespace s t d ;

## typedef SparseGeMatrix<CRS<TYPE> > SPMatrix ;

typedef DenseVector<Array<TYPE> > DEVector ;

SPMatrix ∗W;
DEVector ∗ in ;
DEVector ∗ back ;
DEVector ∗ out ;
DEVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new SPMatrix (N,N ) ;
i n = new DEVector (N ) ;
back = new DEVector (N ) ;
out = new DEVector ( 2 ∗N ) ;
x = new DEVector (N ) ;

## f o r ( i n t i =1; i<=N; ++i )

{
i f ( r a n d v a l () < i c o n n ) ( ∗ i n ) ( i ) = r a n d v a l ( ) ;
i f ( r a n d v a l () < f b c o n n ) ( ∗ back ) ( i ) = r a n d v a l ( ) ;

19
(∗ x ) ( i ) = randval ( ) ;

## f o r ( i n t j =1; j<=N; ++j )

{
i f ( r a n d v a l () < conn ) ( ∗W) ( i , j ) = r a n d v a l ( ) ;
}
}

## f o r ( i n t i =1; i <=2∗N; ++i )

( ∗ out ) ( i ) = r a n d v a l ( ) ;

// f i n a l i z e s p a r s e m a t r i x
W−> f i n a l i z e ( ) ;
}

TYPE compute ( i n t i t e r )
{
// temporary o b j e c t s
DEVector t ;
DEVector : : View out1 = ( ∗ out ) ( ( 1 ,N) ) , out2 = ( ∗ out ) ( (N+1 ,2∗N ) ) ;

while ( i t e r −−)
{
t = ∗x ; // temporary o b j e c t needed f o r BLAS
∗x = a l p h a ∗ ∗ i n + ∗W ∗ t + b e t a ∗ ∗ back ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =1; i<=N; ++i )
( ∗ x ) ( i ) = tanh ( ( ∗ x ) ( i ) ) ;

// c a l c u l a t e y = [ o u t ] ∗ [ x ; x ˆ 2 ]
y = dot ( out1 , ∗x ) ;
f o r ( i n t i =1; i<=N; ++i )
y += out2 ( i ) ∗ pow ( ( ∗ x ) ( i ) , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
DEVector : : View t 1 = ( ∗ x ) ( ( 1 ,N ) ) ;

while ( i t e r −−)
{
∗x = ∗W ∗ t 1 ;
}
}

## Listing 6: Gmm++ Non-Sparse Implementation

// GMM++ NONSPARSE IMPLEMENTATION

#include ”gmm/gmm. h”

#include ” . . / bench−t o o l s . h”

## typedef gmm : : row matrix< s t d : : v e c t o r <TYPE> > DMatrix ;

20
typedef s t d : : v e c t o r <TYPE> DVector ;

DMatrix ∗W;
DVector ∗ in ;
DVector ∗ back ;
DVector ∗ out ;
DVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new DMatrix (N,N ) ;
i n = new DVector (N ) ;
back = new DVector (N ) ;
out = new DVector ( 2 ∗N ) ;
x = new DVector (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
( ∗ i n ) [ i ] = ( r a n d v a l () < i c o n n ) ? r a n d v a l ( ) : 0 ;
( ∗ back ) [ i ] = ( r a n d v a l () < f b c o n n ) ? r a n d v a l ( ) : 0 ;
(∗ x ) [ i ] = randval ( ) ;

## f o r ( i n t j =0; j <N; ++j )

{
( ∗W) ( i , j ) = ( r a n d v a l () < conn ) ? r a n d v a l ( ) : 0 ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) [ i ] = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )
{
// temporary o b j e c t s
DVector t 1 (N ) ;
DVector t 2 ( 2 ∗N ) ;

while ( i t e r −−)
{
gmm : : copy ( ∗ x , t 1 ) ;
gmm : : mult ( ∗W, t1 , ∗x ) ;
gmm : : add ( ∗ x , gmm : : s c a l e d ( ∗ in , a l p h a ) , ∗x ) ;
gmm : : add ( ∗ x , gmm : : s c a l e d ( ∗ back , b e t a ) , ∗x ) ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =0; i <N; ++i )
( ∗ x ) [ i ] = tanh ( ( ∗ x ) [ i ] ) ;

// c a l c u l a t e y = [ o u t ] ∗ [ x ; x ˆ 2 ]
y = gmm : : v e c t s p (gmm : : s u b v e c t o r ( ∗ out , gmm : : s u b i n t e r v a l ( 0 , N) ) , ∗x ) ;
f o r ( i n t i =0; i <N; ++i )

21
y += ( ∗ out ) [ N+i ] ∗ pow ( ( ∗ x ) [ i ] , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
DVector t 1 (N ) ;

while ( i t e r −−)
{
gmm : : copy ( ∗ x , t 1 ) ;
gmm : : mult ( ∗W, t1 , ∗x ) ;
}
}

## Listing 7: Gmm++ Sparse Implementation

// GMM++ SPARSE IMPLEMENTATION

#include ”gmm/gmm. h”

#include ” . . / bench−t o o l s . h”

## typedef gmm : : c s r m a t r i x <TYPE,1 > SPMatrix ;

typedef gmm : : r s v e c t o r <TYPE> SPVector ;
typedef s t d : : v e c t o r <TYPE> DVector ;

SPMatrix ∗W;
SPVector ∗ i n ;
SPVector ∗ back ;
DVector ∗ out ;
DVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new SPMatrix (N,N ) ;
i n = new SPVector (N ) ;
back = new SPVector (N ) ;
out = new DVector ( 2 ∗N ) ;
x = new DVector (N ) ;

// we need temporary w r i t e a b l e s p a r s e o b j e c t s :
gmm : : row matrix< gmm : : w s v e c t o r <double> > W1(N,N ) ;
gmm : : w s v e c t o r <double> i n 1 (N ) ;
gmm : : w s v e c t o r <double> back1 (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
i f ( r a n d v a l () < i c o n n ) i n 1 [ i ] = r a n d v a l ( ) ;
i f ( r a n d v a l () < f b c o n n ) back1 [ i ] = r a n d v a l ( ) ;
(∗ x ) [ i ] = randval ( ) ;

22
f o r ( i n t j =0; j <N; ++j )
{
i f ( r a n d v a l () < conn ) W1( i , j ) = r a n d v a l ( ) ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) [ i ] = r a n d v a l ( ) ;

// now copy t o t h e o r i g i n a l s p a r s e m a t r i x
gmm : : copy (W1, ∗W) ;
gmm : : copy ( in1 , ∗ i n ) ;
gmm : : copy ( back1 , ∗ back ) ;
}

TYPE compute ( i n t i t e r )
{
// temporary o b j e c t s
DVector t 1 (N ) ;
DVector t 2 ( 2 ∗N ) ;

while ( i t e r −−)
{
gmm : : copy ( ∗ x , t 1 ) ;
gmm : : mult ( ∗W, t1 , ∗x ) ;
gmm : : add ( ∗ x , gmm : : s c a l e d ( ∗ in , a l p h a ) , ∗x ) ;
gmm : : add ( ∗ x , gmm : : s c a l e d ( ∗ back , b e t a ) , ∗x ) ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =0; i <N; ++i )
( ∗ x ) [ i ] = tanh ( ( ∗ x ) [ i ] ) ;

// c a l c u l a t e y = [ o u t ] ∗ [ x ; x ˆ 2 ]
y = gmm : : v e c t s p (gmm : : s u b v e c t o r ( ∗ out , gmm : : s u b i n t e r v a l ( 0 , N) ) , ∗x ) ;
f o r ( i n t i =0; i <N; ++i )
y += ( ∗ out ) [ N+i ] ∗ pow ( ( ∗ x ) [ i ] , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
DVector t 1 (N ) ;

while ( i t e r −−)
{
gmm : : copy ( ∗ x , t 1 ) ;
gmm : : mult ( ∗W, t1 , ∗x ) ;
}
}

## Listing 8: MTL Non-Sparse Implementation

// MTL NONSPARSE IMPLEMENTATION

## #include ” mtl / matrix . h”

23
#include <mtl / mtl . h>
#include <mtl / u t i l s . h>

#include ” . . / bench−t o o l s . h”

## typedef matrix<TYPE, r e c t a n g l e <>, dense <>,

row major > : : type DMatrix ;
typedef dense1D<TYPE> DVector ;

DMatrix ∗W;
DVector ∗ in ;
DVector ∗ back ;
DVector ∗ out ;
DVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new DMatrix (N,N ) ;
i n = new DVector (N ) ;
back = new DVector (N ) ;
out = new DVector ( 2 ∗N ) ;
x = new DVector (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
( ∗ i n ) [ i ] = ( r a n d v a l () < i c o n n ) ? r a n d v a l ( ) : 0 ;
( ∗ back ) [ i ] = ( r a n d v a l () < f b c o n n ) ? r a n d v a l ( ) : 0 ;
(∗ x ) [ i ] = randval ( ) ;

## f o r ( i n t j =0; j <N; ++j )

{
( ∗W) ( i , j ) = ( r a n d v a l () < conn ) ? r a n d v a l ( ) : 0 ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) [ i ] = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )
{
// temporary o b j e c t
DVector t 1 (N ) ;

while ( i t e r −−)
{
copy ( ∗ x , t 1 ) ;
mult ( ∗W, t1 , ∗x ) ;
add ( ∗ x , s c a l e d ( ∗ in , a l p h a ) , ∗x ) ;
add ( ∗ x , s c a l e d ( ∗ back , b e t a ) , ∗x ) ;

24
// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =0; i <N; ++i )
( ∗ x ) [ i ] = tanh ( ( ∗ x ) [ i ] ) ;

## y = dot ( ( ∗ out ) ( 0 ,N) , ∗x ) ;

f o r ( i n t i =0; i <N; ++i )
y += ( ∗ out ) [ N+i ] ∗ pow ( ( ∗ x ) [ i ] , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
DVector t 1 (N ) ;

while ( i t e r −−)
{
copy ( ∗ x , t 1 ) ;
mult ( ∗W, t1 , ∗x ) ;
}
}

## Listing 9: MTL Sparse Implementation

// MTL SPARSE IMPLEMENTATION

## #include ” mtl / matrix . h”

#include <mtl / mtl . h>
#include <mtl / u t i l s . h>

#include ” . . / bench−t o o l s . h”

## typedef matrix<TYPE, r e c t a n g l e <>,

compressed <>, row major > : : type SPMatrix ;
typedef dense1D<TYPE> DVector ;
typedef compressed1D<TYPE> SPVector ;

SPMatrix ∗W;
DVector ∗ i n ;
DVector ∗ back ;
DVector ∗ out ;
DVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new SPMatrix (N,N ) ;
i n = new DVector (N ) ;
back = new DVector (N ) ;
out = new DVector ( 2 ∗N ) ;

25
x = new DVector (N ) ;

## f o r ( i n t i =0; i <N; ++i )

{
i f ( r a n d v a l () < i c o n n ) ( ∗ i n ) [ i ] = r a n d v a l ( ) ;
i f ( r a n d v a l () < f b c o n n ) ( ∗ back ) [ i ] = r a n d v a l ( ) ;
(∗ x ) [ i ] = randval ( ) ;

## f o r ( i n t j =0; j <N; ++j )

{
i f ( r a n d v a l () < conn ) ( ∗W) ( i , j ) = r a n d v a l ( ) ;
}
}

## f o r ( i n t i =0; i <2∗N; ++i )

( ∗ out ) [ i ] = r a n d v a l ( ) ;
}

TYPE compute ( i n t i t e r )
{
// temporary o b j e c t
DVector t 1 (N ) ;

while ( i t e r −−)
{
copy ( ∗ x , t 1 ) ;
mult ( ∗W, t1 , ∗x ) ;
add ( ∗ x , s c a l e d ( ∗ in , a l p h a ) , ∗x ) ;
add ( ∗ x , s c a l e d ( ∗ back , b e t a ) , ∗x ) ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =0; i <N; ++i )
( ∗ x ) [ i ] = tanh ( ( ∗ x ) [ i ] ) ;

## y = dot ( ( ∗ out ) ( 0 ,N) , ∗x ) ;

f o r ( i n t i =0; i <N; ++i )
y += ( ∗ out ) [ N+i ] ∗ pow ( ( ∗ x ) [ i ] , 2 ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
DVector t 1 (N ) ;

while ( i t e r −−)
{
copy ( ∗ x , t 1 ) ;
mult ( ∗W, t1 , ∗x ) ;
}
}

## Listing 10: Newmat Non-Sparse Implementation

// NEWMAT NONSPARSE IMPLEMENTATION

## #define WANT STREAM

26
#define WANTMATH

#include ”newmatap . h”
#include ” newmatio . h”

## using namespace NEWMAT;

using namespace s t d ;

#include ” . . / bench−t o o l s . h”

Matrix ∗W;
ColumnVector ∗ i n ;
ColumnVector ∗ back ;
RowVector ∗ out ;
ColumnVector ∗x ;

## i n t main ( i n t argc , char ∗ argv [ ] )

{
run benchmarks ( argc , argv ) ;
}

void i n i t m a t r i c e s ( )
{
W = new Matrix (N,N ) ;
i n = new ColumnVector (N ) ;
back = new ColumnVector (N ) ;
out = new RowVector ( 2 ∗N ) ;
x = new ColumnVector (N ) ;

## f o r ( i n t i =1; i<=N; ++i )

{
i f ( r a n d v a l () < i c o n n ) ( ∗ i n ) ( i ) = r a n d v a l ( ) ;
i f ( r a n d v a l () < f b c o n n ) ( ∗ back ) ( i ) = r a n d v a l ( ) ;
(∗ x ) ( i ) = randval ( ) ;

## f o r ( i n t j =1; j<=N; ++j )

{
i f ( r a n d v a l () < conn ) ( ∗W) ( i , j ) = r a n d v a l ( ) ;
}
}

## f o r ( i n t i =1; i <=2∗N; ++i )

( ∗ out ) ( i ) = r a n d v a l ( ) ;
}

double compute ( i n t i t e r )
{
ColumnVector t ( 2 ∗N ) ;

while ( i t e r −−)
{
∗x = ∗W ∗ ∗x ;
∗x += a l p h a ∗ ∗ i n ;
∗x += b e t a ∗ ∗ back ;

// a c t i v a t i o n f u n c t i o n
f o r ( i n t i =1; i<=N; ++i )

27
( ∗ x ) ( i ) = tanh ( ( ∗ x ) ( i ) ) ;

## // make a temporary v e c t o r ( x , x ˆ2)

t = ∗x & ∗x ;
f o r ( i n t i =1; i<=N; ++i )
t (N+i ) = pow ( ( ∗ x ) ( i ) , 2 ) ;

// compute o u t p u t
y = DotProduct ( ∗ out , t ) ;
}

return y ;
}

## void mt x vec mult ( i n t i t e r )

{
while ( i t e r −−)
{
// STRANGE: t h i s d o e s n o t work :
// x = W ∗ x ;
// i t t a k e s r e a l l y a l o t cpu
// ( b u t i n compute ( ) t h i s works f a s t ! ? )
// so f o r t h e t e s t now t h i s :
∗ i n = ∗W ∗ ∗x ;
}
}

References
 D. Kincaid F. T. Krogh C. L. Lawson, R. J. Hanson. Basic linear subprograms for fortran
usage. 1979. ACM Transactions on Mathematical Software.

##  cplusplus.com. clock reference. http://www.cplusplus.com/reference/clibrary/

ctime/clock.html; accessed 07-08-2007.

##  C.Bischof S.Blackford J.Demmel J.Dongarra J.Croz A.Greenbaum S.Hammarling

A.McKenney D.Sorensen E.Anderson, Z.Bai. Lapack users’ guide. 1999. Society for Indus-
trial and Applied Mathematics.

 Michael Lehn. Everything you always wanted to know about flens, but were afraid to ask.
2007. University of Ulm, Department of Numerical Analysis, Germany.

 Ulisses Mello and Ildar Khabibrakhmanov. On the reusability and numeric efficiency of
c++ packages in scientific computing. 2003. IBM T.J. Watson Research Center, Yorktown,
NY, USA.

##  NC State University. Serial profiling and timing. http://www.ncsu.edu/itd/hpc/

Documents/sprofile.php; accessed 07-08-2007.

28
{"ft_lang_label":"__label__en","ft_lang_prob":0.71765125,"math_prob":0.99194527,"size":18703,"snap":"2019-51-2020-05","text_gpt3_token_len":6290,"char_repetition_ratio":0.14695974,"word_repetition_ratio":0.4902784,"special_character_ratio":0.36502165,"punctuation_ratio":0.17356321,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.99941814,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2019-12-13T02:55:12Z\",\"WARC-Record-ID\":\"<urn:uuid:86174ebd-5335-449b-941e-9d7f287a5ff5>\",\"Content-Length\":\"361873\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:6f732f3d-c286-4092-8afc-f4da0f83d6fb>\",\"WARC-Concurrent-To\":\"<urn:uuid:987d74a0-3f88-4080-bddf-c5774bd1e9d2>\",\"WARC-IP-Address\":\"151.101.202.152\",\"WARC-Target-URI\":\"https://www.scribd.com/document/35851464/SparseLibraryBenchmark-0\",\"WARC-Payload-Digest\":\"sha1:RCGY5W2FCVBVVO353FKHMJXE27B23DIB\",\"WARC-Block-Digest\":\"sha1:DZXZO7VWEHQ3TDZMERE6KE2ERXXFELDS\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2019/CC-MAIN-2019-51/CC-MAIN-2019-51_segments_1575540548537.21_warc_CC-MAIN-20191213020114-20191213044114-00060.warc.gz\"}"}