Global finite element matrix construction based on a
CPU-GPU implementation

arXiv:1501.04784v1 [cs.NA] 20 Jan 2015

Francisco Javier Ramı́rez-Gila , Marcos de Sales Guerra Tsuzukib , Wilfredo
Montealegre-Rubioa,∗
a Department

of Mechanical Engineering, Faculty of Mines, Universidad Nacional de
Colombia, Medellı́n, 050034, Colombia
b Computational Geometry Laboratory, Department of Mechatronics and Mechanical
Systems Engineering, Escola Politécnica da Universidade de São Paulo, Av. Prof. Mello
Moraes, 2231 São Paulo - SP, Brazil

Abstract
The finite element method (FEM) has several computational steps to numerically solve a particular problem, to which many efforts have been directed to
accelerate the solution stage of the linear system of equations. However, the finite element matrix construction, which is also time-consuming for unstructured
meshes, has been less investigated. The generation of the global finite element
matrix is performed in two steps, computing the local matrices by numerical
integration and assembling them into a global system, which has traditionally
been done in serial computing. This work presents a fast technique to construct
the global finite element matrix that arises by solving the Poisson’s equation
in a three-dimensional domain. The proposed methodology consists in computing the numerical integration, due to its intrinsic parallel opportunities, in the
graphics processing unit (GPU) and computing the matrix assembly, due to its
intrinsic serial operations, in the central processing unit (CPU). In the numerical integration, only the lower triangular part of each local stiffness matrix is
computed thanks to its symmetry, which saves GPU memory and computing
time. As a result of symmetry, the global sparse matrix also contains non-zero
∗ Corresponding

author
Email addresses: fjramireg@unal.edu.co (Francisco Javier Ramı́rez-Gil),
mtsuzuki@usp.br (Marcos de Sales Guerra Tsuzuki), wmontealegrer@unal.edu.co (Wilfredo
Montealegre-Rubio )

Preprint submitted to arXiv.org

January 21, 2015

elements only in its lower triangular part, which reduces the assembly operations
and memory usage. This methodology allows generating the global sparse matrix from any unstructured finite element mesh size on GPUs with little memory
capacity, only limited by the CPU memory.
Keywords: Finite element method (FEM), unstructured mesh, matrix
generation, parallel computing, graphic processing unit (GPU), heterogeneous
computing

1. Introduction
Numerous physical phenomena in a stationary condition such as electrical and magnetic potential, heat conduction, fluid flow and elastic problems
in a static condition can be described by elliptic partial differential equations
(EPDEs). The EPDEs does not involve a time variable, and so describes the
steady state response. A linear EPDE has the general form given in Eq. (1),
where a, b, c are coefficient functions, f is a source (excitation) function and u
is the unknown variable. All of these functions can vary spatially (x, y, z).

∇(c · ∇u) + b · ∇u + au = f

(1)

EPDEs can be solved exactly by mathematical procedures such as Fourier
series [1]. However, the classical solution does not frequently exist; for those
problems which allow the use of these analytical methods, many simplifications
are made [2]. Consequently, several numerical methods have been developed
such as the finite element method (FEM) and finite difference method to efficiently solve EPDEs.
The FEM has several advantages over other methods. The main advantage
is that it is particularly useful for problems with complicated geometry using
unstructured meshes [2]. One way to get a suitable framework for solving EPDEs
problems, by using FEM, is formulating them as variational problems, also called
weak solution.

2

The variational formulation of an EPDE is a mathematical treatment for
converting the strong formulation into a weak formulation, which permits the
approximation in elements or subdomains, and the EPDE variational form is
solved using the so-called Galerking method [3]. Fig. 1 summarizes the steps
required to numerically solve an EPDE by using FEM. Some projects are publicly available and all of these steps are solved automatically in the computer
[4]; however, the most common commercial FEM packages are focused on steps
4-8, which are:
(i) Domain discretization with finite elements (FEs);
(ii) Numerical integration of all FEs. In this step, a local matrix ke and a
local load vector f e are computed for each finite element;
(iii) Construction of the global sparse matrix K from local matrices ke and
global load vector F from local load vectors f e ;
(iv) Application of boundary conditions;
(v) Solution of the linear equation system formed previously, KU = F , where
U represents the vector of the nodal solution.

Figure 1: FEM steps to solve EPDEs.

Many efforts have been made to accelerate the linear equation solver by using
direct [5] and iterative methods [6]; however, the FE sparse matrix construction
(steps 5 and 6 of Fig. 1), which are also time-consuming for dense unstructured
meshes, have been less investigated [7]. Hence, in this work the fast construction
3

of the global sparse matrix is developed, which arises in the numerical solution
of EPDEs by FEM.
One way to efficiently construct the global FE matrix is by using parallel
computing. Parallelism is a trend in current computing and, for this reason,
microprocessor manufacturers focus on adding cores rather than on increasing
single-thread performance [8]. Although there are several parallel architectures,
the most widespread many-core architecture are the GPUs (graphics processing units), and the most common multi-core architecture are the CPUs (central
processing units) [9]. However, due to the massive parallelism offered by GPU
compared to CPU, the trend in scientific computing is to use GPU for accelerating intensive and intrinsic parallel computations [10]. Additionally, the GPU
is preferred over the other architectures to be used as a co-processor to accelerate part or the whole code because its market is well established in the
computer games industry. This fact allows having a graphic processor in most
personal computers, which makes GPUs an economic and attractive option for
application developers [10].
Although some researches are focused on completely porting the FEM subroutines to be executed in GPU [11, 12], other works are focused on porting parts
of the FEM code to the GPU, such as numerical integration [13, 14, 15, 16],
matrix assembly [17, 18], and the solution of large linear systems of equations
[19, 20, 21]. Even though some previous works accelerated the construction
of sparse FE matrix [7, 22, 23], these approaches require high GPU memory,
which means that the size of the problem that may be analyzed is limited or an
expensive GPU with high memory capacity is needed. However, in most personal computers equipped with a graphics processor, the GPU has low-memory
capacity.
Accordingly, the proposed approach is based on the assumption that a single GPU has significantly less memory than the CPU; therefore, a CPU-GPU
implementation is developed. The central idea is computing the numerical integration in the GPU, which is an intensive and intrinsic parallel subroutine, and
computing the assembly in the CPU, which is a predominantly serial subroutine
4

with high memory consumption. Hence, based on this methodology, we are able
to construct any large sparse matrix arising in FEM, only limited by the CPU
memory.
The rest of the paper is organized as follows. In section 2, the FEM formulation of the Poisson’s equation used as an example for checking the proposed
implementation is given. Section 3 presents the CPU-GPU implementation of
the FE matrix construction. The results are discussed in section 4 and finally,
in section 5, the conclusions are drawn.

2. Elliptic PDE modelling by FEM
The computational benefits of our proposed implementation in the construction of large sparse FE matrices is shown by solving the Poisson’s equation in a
three-dimensional (3D) domain. The Poisson’s equation is an elliptic-type PDE.
This type of equation is commonly used to solve scalar problems such as electrical and magnetic potential, heat transfer and other problems with no external
flows [3]. In addition, the Poisson’s equation is used for modelling vector field
problems such as structural problems without external applied forces. All of
these problems should be in static or stationary condition. In this section, the
FEM steps are presented, from the strong formulation step to the construction
of the global sparse matrix step (steps 1 to 6 of Fig. 1).
The strong formulation of the Poisson equation is given by [3]:
∇ · (c∇φ(r)) = 0 in Ω
φ = φ̂ on Γφ

(2)

qn = 0 on Γq
where φ represents a general scalar variable on the domain Ω, i.e. can be the
electrical potential or temperature, depending on the specific applications. The
position vector is r and the term c can represent the electrical or thermal conductivity for an isotropic material. Term φ̂ is a prescribed value of scalar variable φ
at boundary Γφ (Dirichlet condition) and the normal component of the flow (or

5

flux) qn is stated for simplicity as zero at boundary Γq (Neumann condition) for
the purposes of computational evaluation only. After the integration by parts
of Eq. (2), the weak form is found as:
Z

(∇ν)T (c∇φ)dΩ = 0

(3)

Ω

where ν is an arbitrary function [3]. Based on the weak form of the Poisson equation, the FE formulation can be found. Accordingly, the domain is discretized
into FEs; next, in each FE, the scalar variable is approximated as:

φ ≈ φ̂ =

n
X

Na φ̃a = N φ̃

(4)

a=1

where n is the number of nodes in the element and, according to the Galerkin
principle, ν = Na are the shape functions. Substituting Eq. (4) into Eq. (3)
and evaluating the integrals for all elements, a system of equation is obtained
as:

KΦ = F ;

F =0

(5)

where K is the so-called global sparse matrix or stiffness matrix. Term Φ is
the nodal solution vector and global load vector F is modified by the boundary
conditions, ceasing to be a vector of zeroes. The K matrix is obtained as a sum
of contributions of the matrices obtained in numerical integration for each FE.
This step is known as assembly [3]:

K=

nel
X

ke

(6)

e=1

where nel is the number of FEs and ke is the local matrix of element e (local
stiffness matrix), which is obtained as:
Z
ke =

cB T BdΩ, with B =

Ωe

n
X

∇Na

(7)

a=1

In our approach we consider that coefficient c, as a part of material properties array (M P ), can vary from one element to the other, but it remains
6

constant for a specific element, which means it can be outside the integral in
Eq. (7). The global sparse matrix construction is made in Eq. (6) and Eq. (7).
However, for computational efficiency, Eq. (7) is solved numerically instead of
by exact integration. There are several numerical procedures to solve an integral, but the most common in FEM is the Gauss-Legendre, or simply Gauss
quadrature [3]. For using the Gauss quadrature, a transformation by using
the reference element technique must be applied, in which physical coordinates
(x, y, z) are transformed into natural coordinates (r, s, t). The equivalence between these systems of coordinates in the local matrix calculation (see the Eq.
(7)) is represented by:
Z
ke

=

c
Ωe
Z1

=

B T (x, y, z)B(x, y, z)dxdydz
Z1 Z1

c

(8)
T

B (r, s, t)B(r, s, t)|J (r, s, t)|drdsdt

−1 −1 −1

where J is the Jacobian matrix and term |J | is its determinant. After that, the
Gauss quadrature is applied to the normalize domain (right hand side of Eq.
(8)) as:

ke = c

Mj Mk
Mi X
X
X

B T (ri , sj , tk )B(ri , sj , tk )|J (ri , sj , tk )|Wi Wj Wk

(9)

i=1 j=1 k=1

where Mi , Mj , Mk , ri , sj , tk and Wi , Wj , Wk are the number, points and weights
of the integration points in each direction in a 3D domain. In our implementation, we use hexahedral FEs with eight nodes in the discretization of the domain,
also known as 8-node brick elements which is a low-order tri-linear FE. The 8node brick element was selected based on the idea that the majority of FE codes
are based on the low-order elements and few works have focused on this type of
elements [13]. The points and weights for an exact numerical integration of the
8-node brick element can be consulted in [3].

7

3. Numerical implementation
A typical way to construct the global sparse FE matrix is using serial computing. There are several ways to construct sparse matrices; however, the fastest
method is to create the triplet1 form first, and then, using a sparse function to
compress the triplet sparse format, avoiding some memory challenges [24]. Despite the triplet format being simple to create, it is difficult to use in most sparse
matrix algorithms; a compress sparse format is thus required [5]. A common
compress sparse format is the CSC (compressed sparse column), used by default
in MATLAB [25].
Algorithm 1 shows a typical implementation for the sparse matrix generation
using the MATLAB syntax. This algorithm needs a subroutine that makes the
numerical integration for computing local matrix ke and local load vector f e .
However, this work is focused only on generating matrix K; load vector F is
not considered for parallel computation (this term is not considered because it
is computationally much less demanding [15, 16]). The numerical integration
subroutine changes with the PDE, element type and basis functions; it is a
function of the nodal coordinates (C) and any nodal fields such as forces or
boundary conditions (BCs) and material properties (M P ) of element e [17].
After computing local matrix ke , which is stored in array K all , a mapping from
local to global degrees of freedom (DOFs) is required. The mapping is made
considering connectivity matrix (E), which gives the global node numbers that
compose an element, and the number of DOFs per node (DOFxn). As a result
of the mapping function, the positions in the global sparse matrix of each entry
of ke is obtained, i.e. row (Indrow ) and column (Indcol ) indices, which are
1 In

the triplet sparse format the row and column indices that specify the position of NNZ

(non-zero) entries must be stored. This means that each NNZ entry requires three numbers:
the row and column position using integer numbers and the value of that entry using a double
precision number. In this format, several NNZ entries with the same row/column indices can
be repeated, which makes the difference with the coordinate (COO) sparse matrix format
where entries with same row/column indices are summed.

8

stored in arrays K rows and K cols , respectively. Finally, the global sparse matrix
is obtained in a single call of an assembly function (e.g. the MATLAB sparse
function [5, 24]).
Algorithm 1 Global sparse FE matrix construction in serial computing.
1: Initialize K rows , K cols and K vals to zero;
2:

for each element e in nel do

3:

ke ←− N umericalIntegration(C, M P, BCs);

4:

K vals (e, :) = ke (:);

5:

(Indrow , Indcol ) ←− M apping(E(e), DOF xn);

6:

K rows (e, :) = Indrow (:);

7:

K cols (e, :) = Indcol (:);

8:

end for

9:

K = sparse(K rows , K cols , K vals );
Several methods have been proposed to parallelize the global sparse matrix

construction. In this work, instead of viewing the matrix generation subroutine
as a single block, it is separated in two subroutines: the numerical integration
subroutine and the assembly subroutine, which allows the analysis looking for
parallelism opportunities. The numerical integration subroutine has several
opportunities for parallelism [16]:
1. Computing all Gauss integration points in parallel;
2. Parallelizing each ke component [15, 16];
3. Finally, it is also possible compute ke in parallel for all elements [13, 16].
The loop over integration points (first approach) is a good strategy for higherorder FEs but it is harder to parallelize due to inherent data race2 in updating
entries in ke with contributions from different integration points, and the degree
2 The

race condition or race hazard is the behavior of a system in which the output is

dependent on a sequence of events. It becomes a bug when events do not happen in the
proper order. Race conditions can occur in computer software, especially multithreaded or
distributed programs.

9

of concurrency is generally lower [16]. The second approach is also helpful for
higher-order FEs, but not convenient for low-order elements. The parallelization
of the loop over low-order elements (third approach) is the most natural way
because a thread can be charged entirely of one local matrix ke [16, 17]. Hence,
the third approach is selected in our implementation since the 8-node brick FE
offers low degree of parallelism using the first two approaches.
In the computation of ke there are no dependences and just the information
from the specific 8-node brick element is necessary. As a consequence, numerical
integration for different elements are also perfectly parallelizable [15]. A thread
is used here in to compute each local matrix ke . Then, with a one-dimensional
grid of thread blocks, all FEs in the mesh can be computed.
In the eight-node brick element for scalar problems (1 DOF per node), each
ke is a square dense matrix having 64 (8 × 8) components, but exploiting the
matrix symmetry, only the lower or upper triangular matrix part, including its
diagonal, should be considered. This leads to a significant saving in floating
point operations and memory requirements, since approximately half the local matrix components are needed, specifically, 36 components instead of 64.
Additionally, there are calculations that need to be made only once for the
same FE type, such as the shape function and their derivatives in the physical
coordinates. We use the same FE type to discretize the whole domain.
Algorithm 2 presents the CUDA kernel to compute the local stiffness matrix
ke . A thread is responsible of computing the 36 components of ke and store
them in the GPU global memory as a vector in array K vals . For each integration point, each thread must serially perform several operations, such as the
Jacobian matrix (J ) computation, its determinant (|J |) and its inverse (J −1 ),
in addition to the gradient matrix (B) and the local stiffness matrix (ke ). The
shape functions (N ) and their derivatives with respect to natural coordinates
(dN ) evaluated at each integration point are computed only once and they are
stored in a fast GPU memory. The computation of matrix J requires three
dense matrix-vector multiplications, i.e. the multiplications of dN (3 × 8) by
nodal coordinates in Cartesian coordinates (8 × 3). The computation of |J | is
10

performed with the determinant formula of a 3 × 3 matrix and J −1 with the
formula to invert matrices of the same size. Matrix B is computed as the multiplication of J −1 (3 × 3) by dN (3 × 8). Finally, element matrix ke is computed
as presented in Eq. (9) with some modifications to calculate only the lower
triangular matrix components.
Algorithm 2 CUDA kernel for numerical integration in the 8-node brick element.
1:
global

void N umericalIntegration(ke , E, C, M P )

2:

tid = blockDim.x ∗ blockIdx.x + threadIdx.x; // Thread ID

3:

// Initialize local variables

4:

if tid < nel then

5:

for i = 0; i < 8; i + + do // Loop over integration points

6:

read N (i) and dN(i);

7:

compute J , |J | and J −1 ;

8:

compute B;

9:

compute the lower triangular part of ke and store it in K vals ;

10:
11:

end for
end if
However, if the FE mesh is denser, the numerical integration may not be

executed because this is a heavily memory-consuming operation and the GPU
memory cannot allocate it entirely [17]. Thus, a simple domain partitioning
based on the GPU memory availability is used, which is a modified implementation of that proposed by Komatitsch et al. [11]. Consequently, the numerical
integration is developed in subgroups of FEs (GP Ug ), i.e. the numerical integration kernel is invocated to operate a FEs set. The elements subgroups
are executed sequentially by the GPU, but the numerical integration over FEs
within a group are executed in parallel. With this approach, the construction
of the global FE matrix can be made from any mesh size, only limited by the
CPU memory. Thus, if the available GPU memory is GP Uma and the required
memory by the numerical integration subroutine is GP Umr (known through
11

GPU query functions), the work done by the GPU is divided into groups based
on the following formula:

GP Ug = ceil

GP Umr
GP Uma


(10)

Each group of FEs is integrated numerically in the GPU, and the results
allocated in the global GPU memory are downloaded to the CPU memory,
leaving the GPU memory free for a new group of FEs. Every kernel is uploaded
with only the data required by the FEs belonging the group that will be executed
(only part of arrays C and E), which allows a better use of memory.
After the numerical integration over all elements, the assembly of all local
matrices into the global matrix is executed. To do this, the assembly step
can be viewed as a sparse matrix format conversion [7, 24]. As stated before,
the global position of each component of the K all matrix should be stored in
arrays K rows and K cols , which is an implicit sparse matrix format storage
[5]. The assembly step consists in converting the triplet format into the CSC
format. This function requires sorting and reducing operations, which are not
well suited for parallel computing due to they are intrinsic serial operations. In
fact, the assembly cannot be considered an embarrassingly parallel algorithm
with no dependencies, as is considered the numerical integration is [15]. The
dependencies present in the assembly algorithm can cause race condition, which
should be avoided to prevent wrong results.
To remove dependencies present in the assembly algorithm there are some
options, but the best for GPU programming is the use of coloring techniques or
atomic operations [7, 17]. Using the first option, sets of elements with different
colors are obtained, with two elements having the same color if their local stiffness matrices do not coincide at any global matrix entry. Given this element
partition, the numerical integration algorithm followed by immediate matrix assembly can be executed in parallel for elements within a group of the same color
and proceeds color by color sequentially. However, this option requires precomputing operations to divide the mesh into colors, where the problem to find

12

the optimal number of colors is NP-hard. Although there are many heuristics
for finding nearly minimal colorings [26], this operation can take substantial execution time in meshes with a large number of FEs. On the other hand, atomic
operations are undesirable because the penalty on the GPU execution time can
be high [17].
Dziekonski et al. [7] use the sparse matrix conversion concept, which is executed in the GPU with the advantage that the assembly occurs in parallel using
atomic operations. Using this approach, arrays K rows and K cols must be computed and stored in GPU, which substantially increases memory requirements.
Although this approach is good on high-end GPUs provided with high-memory
capacity to avoid memory transactions between processors, it is less interesting
for many low-end video cards users. Additionally, in many current personal
computers, the memory installed in the CPU motherboard is often four to eight
times as large as than that installed on the graphics card [11]. Hence, a good
option is to make the sparse matrix format conversion in the CPU.
Therefore, a methodology for using GPU for intrinsic parallel and CPU for
intrinsic serial algorithms in the global FE matrix construction is proposed.
The components of K vals are parallel computed in GPU, but indices K rows
and K cols are computed and stored in the CPU, since their computation is
simple [24]. After computing the global FE matrix in triplet format in CPU
and GPU, as described above, the sparse matrix format conversion is executed
in MATLAB using one of the following options: the native MATLAB function
sparse 3 , the sparse2 function from SuiteSparse4 and the sparse create function
from MILAMIN5 . A graph showing the flow of our proposed implementation is
presented in Fig. 2.
In Fig. 2, there is a preprocessing phase that can be performed by any
specialized software for complex domains. However, as our goal is evaluating
3 MATLAB

R2013b, http://www.mathworks.com/help/matlab/ref/sparse.html
4.0.2, http://www.cise.ufl.edu/research/sparse/SuiteSparse/
5 Mutils 0.3, http://www.milamin.org/
4 SuiteSparse

13

the computational time involved in the matrix generation, we use cubic domains
discretized by 8-node brick elements. This simple discretization is performed in
MATLAB generating three arrays: connectivity matrix E, which indicates the
nodes that compose an element; nodal matrix coordinates C, which specify the
Cartesian coordinates of the nodes; and material properties matrix M P , which
states the material properties of each element. As this subroutine is simple and
takes a few seconds to generate large meshes, it is executed in the CPU.
The second phase in Fig. 2 shows the necessary steps for building the global
FE matrix, as described previously. In this step, the heterogeneous computing
concept is used: while the GPU is executing the numerical integration kernel,
the CPU is used for generating the indices of the sparse matrix in triplet format.
Once the local matrices from all elements are computed, the assembly function
can be executed, and the global matrix is obtained. After the global system is
obtained, the boundary conditions are applied, and the solution of the linear
system is found, but these steps are not shown in the figure because it is out of
the scope of this paper.

Figure 2: Global sparse FE matrix construction.

4. Results
The GPU results are obtained using a 2GB NVIDIA Quadro 4000 video
card based on Fermi architecture. The CPU computations were performed in
14

MATLAB R2013b using an Intel Xeon E5-2620 CPU at 2.0 GHz with six cores
in a 64-bit Windows 7 operating system with 40GB of RAM. However, in the
serial implementation only one CPU core is used. As stated before, the matrix
is generated in two steps: first, conducting the numerical integration (NI) over
all FEs, and second, assembling all local matrices into a global matrix. The
results are presented separately in the following subsections.
4.1. Numerical integration
The NI subroutine is programmed in CUDA C and is called from MATLAB workspace using the parallel.gpu.CUDAKernel constructor, which creates
a CUDAKernel object. The properties of the CUDAKernel object such as the
number of blocks in the grid and the number of threads in the block are also
configured in MATLAB. After the creation and configuration of the kernel, it
is evaluated using the feval MATLAB function.
The code executed in GPU is compared with its respective serial implementation counterpart executed in one CPU core. Fig. 3 shows the computational
time when the NI routine is executed in CPU and GPU for different meshes. In
the time spent by the GPU, the time spent in memory transactions is also taking
into account: uploading data from CPU to GPU memory using the MATLAB
function gpuArray, and vice versa, downloading data from GPU to CPU memory using the function gather. The gather function guarantees that all work is
finished in the GPU, and it ensures accurate time measurement using the tic-toc
MATLAB functions. Fig. 3 also shows the speedup reached by the GPU implementation. The GPU speedup is calculated as the ratio between the serial (ts )
and parallel (tp ) execution time. The advantage of the GPU implementation is
observed in the curves presented.
The NI routine can compute meshes with almost five million of 3D FEs in a
GPU with-low-memory capacity without using mesh partitioning. Although the
NI routine is considered memory bounded [17], the implemented code consumes
little memory due to the low-order of the FE selected. The maximum time
spent in the GPU for the largest mesh is 7.56 seconds, while for the same mesh
15

Figure 3: Execution time in the numerical integration subroutine and GPU speedup as a
function of the mesh size.

in serial computing in the CPU, the NI routine spent 875.1 seconds, which
means a GPU speedup of 115.8X. The GPU speedup dramatically increases
its value using 8000 FEs. After that, the GPU speedup rate is slow; indeed,
with more than four million FEs, the GPU speedup decreases. This fact is
explained by the ratio between floating point operations (flops) and memory
operations (memops). When the data volume transferred from one to another
memory type (memops) is larger than the flops; the problem thus becomes
memory bounded and the speedup is decreased. This is known as memory wall,
in which there is a difference between the speed at which data are transferred to
the processor and the rate at which instructions are executed [27]. A common
problem in transferring data between the CPU and the GPU memories is found
in the PCI-express port, which limits the rate of moving data between these
processors. For this reason, minimizing the transaction between the processor
memories is recommended.

16

4.2. Assembly
The assembly step is performed in the CPU for two reasons: the assembly algorithm presents few parallel opportunities and requires high temporary
memory. As the CPU cores are designed to maximize the execution speed of
sequential programs (latency6 ) [10], this processor is well suited for serial tasks,
as proposed here (assembly). The memory requirements of the assembly routine
are shown in Table 1.
As indicated in Section 3, the assembly phase can be viewed as a sparse
matrix conversion from triplet to CSC format. Table 1 shows the percentage of
compression in the non-zero elements (NNZ) going from one to another sparse
format through different meshes. This table also shows the percentage of memory saving when the CSC is used instead of the triplet form. The percentage
of compression and memory saving are not equal due to the matrix structure,
which depends on the nodes and elements numeration; however, these percentages seem to be very similar.
Having into account that the assembly subroutine initially requires a lot of
memory for storing the three arrays K rows , K cols and K vals , it is not feasible
for a GPU implementation, particularly if GPUs with low-memory capacity are
used for forming large matrices.
Three different routines for changing the sparse matrix format are considered
herein: the sparse function from MATLAB (Sparse), the sparse2 function from
SuiteSparse (Sparse2) and the sparse create function from MILAMIN (SparseCreate). The execution time spent in the NI routine, executed in the GPU,
and the assembly routine, executed in the CPU, are presented in Fig. 4. Addi6 Latency,

in a basic definition, is the amount of time to complete a task, which mean that

it is measured in units of time. While the CPU tries to minimize the latency, the GPU tries to
maximize the throughput. Throughput is tasks completed per unit time measured in units as
stuff per time, such as jobs completed per seconds. Therefore, the NI subroutine is executed in
the GPU because we need to obtain many local matrices in the unit of time; and the assembly
subroutine is executed in the CPU because the whole task is required to be finished in a short
time.

17

Table 1: Memory requirements in the triplet and CSC sparse formats for different mesh
discretizations.

# of FEs

1 000

8 000

64 000

512 000

1 728 000

2 744 000

4 096 000

5 832 000

8 000 000

NNZ in Triplet

36 000

288 000

2 304 000

18 432 000

62 208 000

98 784 000

147 456 000

209 952 000

288 000 000

NNZ in CSC

15 561

118 121

920 241

7 264 481

24 408 721

38 710 841

57 728 961

82 135 081

112 601 201

NNZ compression

56.8%

59.0%

60.1%

60.6%

60.8%

60.8%

60.9%

60.9%

60.9%

Memory used (triplet)

0.58 MB

4.61 MB

36.9 MB

294.9 MB

995.3 MB

1 580.5 MB

2 359.3 MB

3 359.2 MB

4 608.0 MB

Memory used (CSC)

0.26 MB

1.96 MB

15.3 MB

120.5 MB

404.7 MB

641.8 MB

957.0 MB

1 361.6 MB

1 866.6 MB

54.9%

57.4%

58.6%

59.1%

59.3%

59.4%

59.4%

59.5%

59.5%

Memory saving

tionally, this figure shows the time spent for computing arrays K rows and K cols
(Index) used by Sparse and Sparse2 routines. The NI routine spends more time
than Sparse2 and SparseCreate routines for all mesh discretizations, even though
the NI routine is executed in parallel in the GPU and the assembly routines are
executed serially in the CPU. This occurs because the NI is an intensive routine
that requires many operations between small matrices, such as multiplications,
additions and inversions, whereas assembly only requires sorting and reducing
operations [7].

Figure 4: Execution time spent by the assembly using three functions and different mesh sizes.

The sparse MATLAB function consumes significantly more time and mem18

ory than the other assembly functions, and for meshes larger than 200 thousand
NNZs, it requires more time than the NI routine. The large memory usage
explains the drawbacks of this function. In this routine, the row and column
indices must be entered as double type arrays, which increases the memory requirements and performance overhead making this function the most inefficient
assembly routine [24]. In turn, although the sparse2 function has an excellent
performance, this routine requires computing and storing the index arrays, and
their computations require almost the same time as the assembly function itself,
meaning that the computational time spent by sparse2 can be duplicated if the
time for computing K rows and K cols is considered.
Regarding the sparse create function for the sparse matrix conversion, in
all discretization cases, the computational time invested in assembly is always
the lowest. Indeed, this function employs less time than the index arrays efficiently computed as MILAMIN [24]. Moreover, the sparse create function does
not require the computation of the index arrays for the assembly step, saving
memory and computations. However, this function is a specialized routine to
form only matrices arising in FEM, which can represent a disadvantage when
general sparse matrices are required.
4.3. Partitioning the numerical integration
Finally, as shown above, the numerical integration executed in parallel by using GPU and assembling the NNZ element based on sparse create CPU function
are the most efficient implementation. However, without doing a mesh partitioning for the NI routine, the problem size that can be solved using low-memory
capacity GPUs is limited. Hence, using an approach as the one presented in
Section 3 for partitioning the mesh based on available GPU memory, FE matrices can be constructed as large as CPU memory allows. In addition, as the
CPU memory is usually greater than the GPU memory, larger matrices can be
constructed than those created only in the GPU.
Fig. 5 shows the computational time for numerical integration and assembly
routines in the construction of the global sparse matrix, considering different
19

mesh sizes. Adding the time spent by these two routines, the time necessary to
construct the global sparse matrix (MatGen) is obtained. Table 2 presents these
results for different mesh sizes and the time percentage employed per function.
The largest mesh size constructed has 27.27 million nodes and 27 million 3D
elements. As stated before, the numerical integration routine always consumes
more time than the assembly. Specifically, the assembly function consumes less
than the 40% of the total matrix generation time for all the mesh cases: 37.9% is
the greatest percentage spent for assembling 1 728 million FEs (an intermediate
mesh size) and the lowest percentage was presented assembling 1 000 FEs (the
lowest mesh size). These results indicate an appropriate approach to generate
large FE matrices from unstructured meshes in systems in which the GPU
memory is lower than the CPU memory. Additionally, this method shows the
benefits of heterogeneous computing: using the CPU for intrinsic serial routines
and the GPU for intrinsic parallel and intensive computations.

Figure 5: Execution time spent generating global sparse FE matrices of several sizes.

As a final comment, for reducing the time for matrix generation, one option
is to form temporally small sparse matrices and next, to assemble all these
matrices for forming the global matrix. Since the numerical integration for large

20

Table 2: Computational time spent in the sparse matrix construction.

# of FEs

1 000

8 000

64 000

512 000

1 728 000

2 744 000

4 096 000

5 832 000

8 000 000

15 625 000

27 000 000

Matrix size

1 331

9 261

68 921

531 441

1 771 561

2 803 221

4 173 281

5 929 741

8 120 601

15 813 251

27 270 901

MatGen time (s)

0.02

0.03

0.15

1.17

3.98

5.92

8.85

15.3

23.8

40.9

115.0

NI time

82.8%

80.9%

68.4%

62.9%

62.1%

66.2%

66.2%

70.9%

67.8%

70.5%

82.4%

Assembly time

17.2%

19.1%

31.6%

37.1%

37.9%

33.8%

33.8%

29.1%

32.2%

29.5%

17.6%

meshes is divided into groups, which are executed sequentially in the GPU, it
can can execute kernels asynchronously7 ; an approach that can accelerate large
sparse matrix generation in FEM is assembling all local matrices belonging to
a FEs group while the GPU is executing the NI routine.

5. Conclusions
The implementation proposed shows that large global sparse matrices arising
in the numerical solution of EPDEs by FEM can be generated in heterogeneous
computing with limited computational resources. In fact, with only 2GB of
memory in the GPU, FEM matrices of 3D unstructured meshes with up to
27 million of FEs are generated. Moreover, a significant acceleration in the
numerical integration executed in the GPU as compared to the serial CPU implementation is achieved. The maximum speedup reached by the GPU is 126x,
and the mean value is 121x, which shows good performance of the developed
code. Additionally, the assembly function (sparse create) used in this work is
an efficient routine that requires lower runtime and memory consumption compared with the other two functions (sparse and sparse2 ). In fact, sparse create
consumes less runtime than the numerical integration routine, even though the
latter is implemented in parallel computing. The maximum time spent by the
sparse create function took only twenty seconds for assembling 27 million FEs,
7 Some

function calls are asynchronous to facilitate concurrent execution, returning the

control to the CPU thread before the GPU has completed a specific task. This means that
the GPU is able to execute several kernels while the CPU executes other routines.

21

one-fifth of the time spent by NI routine.
Finally, although the proposed methodology is applied to the numerical solution of the Poisson equation by using FEM, this approach is general and can
be used to construct not only the stiffness matrix from static or stationary problems, but also the mass matrix from dynamics problems and for solving other
types of PDEs by using the FEM.

Acknowledgments
The first author thanks to Universidad Nacional de Colombia for the financial support to this work with the Estudiantes sobresalientes de posgrado scholarship. MSG Tsuzuki was partially supported by CNPq (grant 310.663/2013-0).

References
References
[1] E. Kreyszig, H. Kreyszig, E. Norminton, Advanced engineering mathematics, 10 ed., John Wiley & Sons, Hoboken, N.J., 2011.
[2] D. Braess, Finite elements: Theory, fast solvers, and applications in solid
mechanics, 2 ed., Cambridge University Press, Cambridge, 2001.
[3] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The finite element method: its
basis and fundamentals, 6 ed., Elsevier, Amsterdam, 2005.
[4] A. Logg, K.-A. Mardal, G. Wells (Eds.), Automated solution of differential
equations by the finite element method: The FEniCS Book, volume 84 of
Lecture Notes in Computational Science and Engineering, Springer, 2012.
[5] T. A. Davis, Direct methods for sparse linear systems, Fundamental of
Algorithms, SIAM, Philadelphia, 2006.
[6] Y. Saad, Iterative methods for sparse linear systems, 2 ed., SIAM, 2003.
doi:10.1137/1.9780898718003.

22

[7] A. Dziekonski, P. Sypek, A. Lamecki, M. Mrozowski, Progress In Electromagnetics Research 128 (2012) 249–265. doi:10.2528/PIER12040301.
[8] J. Owens, M. Houston, D. Luebke, S. Green, J. Stone, J. Phillips, Proceedings of the IEEE 96 (2008) 879–899. doi:10.1109/JPROC.2008.917757.
[9] A. R. Brodtkorb, C. Dyken, T. R. Hagen, J. M. Hjelmervik, O. O. Storaasli,
Scientific Programming 18 (2010) 1–33. doi:10.3233/SPR-2009-0296.
[10] D. B. Kirk, W.-m. W. Hwu, Programming massively parallel processors:
A hands-on approach, 2 ed., Morgan Kaufmann, Burlington, MA, USA,
2012.
[11] D. Komatitsch, D. Michéa, G. Erlebacher, Journal of Parallel and Distributed Computing 69 (2009) 451–460. doi:10.1016/j.jpdc.2009.01.
006.
[12] Z. Fu, T. J. Lewis, R. M. Kirby, R. T. Whitaker, Journal of computational
and applied mathematics 257 (2014) 195–211. doi:10.1016/j.cam.2013.
09.001.
[13] M.

G.

ematical

Knepley,

A.

R.

Terrel,

Software

39

(2013)

1–13.

ACM

Transactions

on

Math-

doi:10.1145/2427023.2427027.

arXiv:arXiv:1103.0066v1.
[14] M. Pawel, P. Plaszewski, K. Banaś, Procedia Computer Science 1 (2010)
1093–1100. doi:10.1016/j.procs.2010.04.121, 00014.
[15] F. Krużel, K. Banaś, Computers & Mathematics with Applications 66
(2013) 2030–2044. doi:10.1016/j.camwa.2013.08.026.
[16] K. Banaś, P. Plaszewski, P. Maciol, Computers & Mathematics with Applications 67 (2014) 1319–1344. doi:10.1016/j.camwa.2014.01.021.
[17] C. Cecka, A. J. Lew, E. Darve, International Journal for Numerical Methods in Engineering 85 (2011) 640–669. doi:10.1002/nme.2989.

23

[18] G. R. Markall, A. Slemmer, D. A. Ham, P. H. J. Kelly, C. D. Cantwell,
S. J. Sherwin, International Journal for Numerical Methods in Fluids 71
(2013) 80–97. doi:10.1002/fld.3648.
[19] J. Krüger, R. Westermann, in: M. Pharr, R. Fernando (Eds.), GPU gems
2: programming techniques for high-performance graphics and generalpurpose computation, number 2 in GPU gems, Addison-Wesley, Upper
Saddle River, NJ, 2005, pp. 703–718. 00467.
[20] F. Stock, A. Koch, in: R. Wyrzykowski, J. Dongarra, K. Karczewski,
J. Wasniewski (Eds.), Parallel Processing and Applied Mathematics, number 6067 in Lecture Notes in Computer Science, Springer Berlin Heidelberg,
2010, pp. 457–466. 00002.
[21] R. Li, Y. Saad, The Journal of Supercomputing 63 (2012) 443–466. doi:10.
1007/s11227-012-0825-3.
[22] A. Dziekonski, P. Sypek, A. Lamecki, M. Mrozowski, International Journal
for Numerical Methods in Engineering 94 (2013) 204–220. doi:10.1002/
nme.4452.
[23] A. Dziekonski, P. Sypek, A. Lamecki, M. Mrozowski, IEEE Antennas
and Wireless Propagation Letters 11 (2012) 1346–1349. doi:10.1109/LAWP.
2012.2227449.
[24] M. Dabrowski, M. Krotkiewski, D. W. Schmid, Geochemistry, Geophysics,
Geosystems 9 (2008) 1–24. doi:10.1029/2007GC001719.
[25] J. R. Gilbert, C. Moler, R. Schreiber, SIAM Journal on Matrix Analysis
and Applications 13 (1992) 333–356. doi:10.1137/0613024.
[26] M. Kubale, Graph Colorings, American Mathematical Society, Providence,
Rhode Island, 2004.
[27] K. Asanovic, J. Wawrzynek, D. Wessel, K. Yelick, R. Bodik, J. Demmel,
T. Keaveny, K. Keutzer, J. Kubiatowicz, N. Morgan, D. Patterson, K. Sen,
Communications of the ACM 52 (2009) 56. doi:10.1145/1562764.1562783.
24

