# Evaluating Accuracy and Efficiency of HPC Solvers for Sparse Linear Systems
with Applications to PDEs

Antonella Galizia⋆ Simone Cammarasana⋆ Andrea Clematis  Giuseppe Patané

###### Abstract

Partial Differential Equations (PDEs) describe several problems relevant to
many fields of applied sciences, and their discrete counterparts typically
involve the solution of sparse linear systems. In this context, we focus on
the analysis of the computational aspects related to the solution of large and
sparse linear systems with HPC solvers, by considering the performances of
direct and iterative solvers in terms of computational efficiency,
scalability, and numerical accuracy. Our aim is to identify the main criteria
to support application-domain specialists in the selection of the most
suitable solvers, according to the application requirements and available
resources. To this end, we discuss how the numerical solver is affected by the
regular/irregular discretisation of the input domain, the discretisation of
the input PDE with piecewise linear or polynomial basis functions, which
generally result in a higher/lower sparsity of the coefficient matrix, and the
choice of different initial conditions, which are associated with linear
systems with multiple right-hand side terms. Finally, our analysis is
independent of the characteristics of the underlying computational
architectures, and provides a methodological approach that can be applied to
different classes of PDEs or with approximation problems.

Keywords. HPC, Numerical Solvers, Sparse linear systems, SuperLU, PETSc, PDE

##  1 Introduction

Partial differential equations (PDEs) are crucial for the understanding of the
behaviour of several phenomena [18] in engineering [14], geophysical
exploration [24], and fluido-dynamics [38]. Since it is not always possible to
compute the analytic solution of PDEs, numerical approaches discretise PDEs
and generally require the solution to linear systems [33]. In this context,
previous work typically applies a workflow that includes the following steps:
modelling of the input problem, discretisation of the input domain and of the
PDE, and solution of the discrete equation, which typically results in a
large, sparse linear system (Fig. 1).

With respect to this framework, we focus on the analysis of computational
aspects related to the solution of large and sparse linear systems with HPC
solvers. We analyse the performance of direct and iterative solvers by taking
into account the impact of regular/irregular discretisation grids of the input
domain on the computational cost and on the numerical accuracy of the
solution. As additional elements of our analysis, we discuss the influence of
the discretisation of the input PDE with piecewise linear or polynomial basis
functions and with different initial conditions. In fact, the discretisation
of PDEs with Finite Element Methods (FEMs) and polynomial basis functions
generally results in a reduced sparsity of the coefficient matrix with respect
to the linear basis, and different initial conditions are associated with
linear systems with multiple right-hand side terms. As test cases, we consider
the Laplace equation on a 2D/3D domain, as it is a standard problem for
testing numerical solvers, and it is related to several physical applications
(e.g., electrostatic potential, heat diffusion in steady state). However, a
similar analysis can be extended to different classes of PDEs, such as
Laplacian eigen problems, heat equation, etc.

To solve sparse linear systems associated with PDEs, we employ parallel
software freely accessible (SuperLU, PETSc [27, 9]) and consider several
configurations, performance metrics, and properties, such as efficiency,
scalability, and accuracy. Instead of evaluating the impact of the underlying
architecture on the performance of the solvers, our objective is to identify
the main criteria to support application specialists in the selection of the
most suitable solvers, according to the underlying application and available
computational resources. Following the approach proposed in [7], our aim is to
derive a “rule of thumbs” for _production-oriented programmers_ , i.e., domain
specialists exploiting, at their best, the tools provided by _optimisation-
oriented programmers_ , who are committed towards the design of efficient and
optimised libraries. With this goal in mind, we characterise sparse linear
solvers by taking into account application constraints, such as the
discretisation of the input domain (e.g., regular/irregular), the properties
of the coefficient matrix (e.g., size, sparsity, conditioning) and of the
solver (e.g., accuracy, scalability). The proposed analysis will support the
exploitation of parallel tools, as available off-the-self solvers, and
possibly reduce the barrier to entry in HPC ecosystems. Furthermore, thanks to
its generality, our discussion represents a methodological approach to be
applied with different classes of PDEs (e.g., elliptic, parabolic PDEs) or
with approximation schemes.

#### Related work, novelty, and main contributions

The scientific community has a great interest in the numerical solution of
linear systems, as confirmed by a long-standing tradition in software
development [17], [10], [36]. An overview on the state-of-the-art of
widespread tools is presented in Sect. 2.1.

(a) | (b) | (c)
---|---|---

(d) | (e)
---|---

Figure 1: (a) Input CAD model, (b) color-map and (c) values of the solution to
the Laplace equation. (d,e) Solution on two orthogonal planes in the inner
volume.

Different works propose efficient parallel solvers of sparse linear systems:
novel strategies and related performance figures are compared with well-known
software packages, together with the discussion of theoretical complexity and
accuracy, and/or the use of specific parallel paradigms, run-time systems and
accelerators. In [13], a parallel hierarchical solver is proposed and compared
with the SuperLU performances. In [39], the performances of a novel multi-
frontal solver are discussed and compared with MUMPS [5] and SuperLU. In [2,
3], a DAG-based parallelisation of the QR factorization of sparse matrices is
implemented over the StarPU run-time system. The GHOST toolkit [25] implements
a collection of building blocks for sparse linear algebra on heterogeneous
systems, such as multicore processors, graphics processing units (GPUs), and
other accelerators (e.g., Intel Xeon Phi). Focusing on GPUs, in [30], NVIDIA
GPU-accelerated libraries are used to implement the incomplete-LU and Cholesky
preconditioned iterative solvers; in [1] and [28], optimised GPU
implementations of key operations in sparse linear systems solvers are
proposed. In [16], a new benchmark for the computation of the _High-
Performance Conjugate Gradient_ for the solution of the Poisson equation on a
regular 3D grid is presented. This work also discusses a closer relation of
the benchmark with application requirements to drive future systems’ design.

This previous work emphasises the need of effective sparse linear solvers from
the developer’s viewpoint, i.e., optimisation-oriented programmers
investigating the development of optimised solvers with notable performances,
the exploitation of accelerators, and the design of future systems. From a
different perspective, our aim is to address the production-oriented
programmers needs through a comparison of well-known tools and the
identification of the main criteria for the selection of the direct or
iterative solvers of sparse linear systems that better suit application
requirements.

More close to our approach, in [32, 12] several parallel solvers of sparse
linear systems discretising PDEs are applied for applications in geophysical
exploration and fluid-dynamics, also discussing solvers’ effectiveness with
respect to the specific applications. In [12], a parallel software environment
for 3D fluid dynamic simulations is presented; three case studies are proposed
to test direct and iterative solvers on different levels of simulation
complexity. While these works discuss solvers for specific problems, we aim at
addressing the needs of a larger class of PDEs. To the best of our knowledge,
an up-to-date comparison of the computational efficiency and numerical
accuracy of direct and iterative solvers of large sparse linear systems (e.g.,
associated with discrete PDEs) is actually missing.

#### Paper organisation

Firstly (Sect. 2), we present the selected case study, the parallel tools
employed to solve the corresponding sparse linear systems, data sets with
different matrix sparsity, metrics and test-beds. Then, we present the
performance figures of direct, and iterative solvers (Sect. 3), which provide
a fine granularity of analysis about scalability, and possible bottlenecks of
each solver. Solvers are then evaluated in terms of memory requirements,
numerical accuracy of the solution, multiple right-hand side vectors and
linear/polynomial basis functions (Sect. 4), also presenting their main
limitations. Finally, (Sect. 5) we compare numerical solvers and we present
the rule of thumbs to support application specialists; last considerations and
future work conclude the proposed analysis (Sect. 6).

##  2 Experimental pipeline

We briefly review the main aspects related to direct and iterative solvers
(Sect. 2.1), our case studies (Sect. 2.2), metrics and testbeds (Sect. 2.3).

###  2.1 Direct and iterative solvers

#### Numerical libraries

Numerical linear algebra has always represented a central topic for the
scientific community; many high-performance software packages and tools for
parallel heterogeneous architectures have been developed in the last decades
[17], [10], [15]. Several studies propose general purpose solvers and
preconditioners of sparse linear systems on parallel resources. For instance,
PaStiX [21], SuperLU [26], MUMPS [5], and PARDISO [35] provide parallel sparse
direct solvers. Hypre [19] is a parallel suite of preconditioners and solvers
for sparse linear systems. Trilinos [22] and PETSc - Portable, Extensible
Toolkit for Scientific Computation [9] represent frameworks for the solution
of complex scientific problems; among others, both libraries provide several
modules and external packages for sparse iterative and direct solvers.

Since all the aforementioned software tools have been extensively used,
updated, and tested on a large set of problems, they represent the ideal base
for our discussion, as we are interested in a systematic analysis of the
performance and accuracy of available parallel solvers for sparse linear
systems. To this end, we experiment different algorithms and tools [11],
compare MUMPS [5, 6] and SuperLU [27] for direct solvers, and focus on PETSc
[9] for iterative solvers and preconditioners.

#### Direct and iterative linear solvers

To discuss performance and accuracy figures, we briefly review the main
aspects and features of direct and iterative solvers of a linear system
$\textbf{A}\textbf{u}=\textbf{b}$. _Direct solvers_ [20] are based on the LU
_decomposition_ , which factorises the coefficient matrix as
$\textbf{A}=\textbf{L}\textbf{U}$, where L and U are a lower and an upper
triangular matrix, respectively. The triangular linear systems
$\textbf{L}\textbf{y}=\textbf{b}$, $\textbf{U}\textbf{u}=\textbf{y}$ are
solved with a forward and a backward substitution in linear time with respect
to the number of unknowns, thus reducing the quadratic or cubic computational
cost for the solution to the input linear system.

If A is sparse, then the factorisation typically introduces many more nonzero
entries, i.e., the $(i,j)$ entry of $\mathbf{A}$ can be non-zero in the
$(\mathbf{L},\mathbf{U})$ factors even if it is originally zero. This
phenomenon is known as _fill-in_ and has an impact on (i) the growth of the
memory requirement and the number of operations, (ii) an a-priori unknown
structure of the factorised matrices. Re-ordering methods (e.g., approximate
minimum degree method [4], multilevel graph partitioning [23]) can be applied
to the columns of the coefficient matrix to minimise the fill-in introduced by
the LU factorisation. Moreover, a symbolic factorisation can be used to define
a sparsity pattern of the LU decomposition with a reduced computational
complexity with respect to the numerical factorisation. To efficiently
generate the (L, U) matrices, we can apply the left-looking (Spooles [8]),
right-looking (SuperLU_DIST), and multi-frontal (MUMPS) method.

_Iterative solvers_ [20, 34] (e.g., Krylov subspace projection, Arnoldi or
Lanczos iterations, transpose-free variants) compute an approximated solution
of a linear system, starting from an initial guess and without factorising the
input coefficient matrix. The efficiency of a solver depends mainly on the
number of operations for each iteration and on the number of iterations needed
to converge to the solution under a certain stopping criterion. The number of
iterations depends also on the condition number of the coefficient matrix,
which can be reduced with a preconditioning operation. We briefly recall that
preconditioners are matrices used to rewrite the input linear system in an
equivalent formulation, whose new coefficient matrix has a lower conditioning.
Given a linear system $\textbf{A}\textbf{u}=\textbf{b}$ and a preconditioner
$\mathbf{M}^{-1}$, we solve the linear systems

$\textbf{A}\mathbf{M}^{-1}\textbf{y}=\textbf{b},\qquad\textbf{u}=\mathbf{M}^{-1}\textbf{y},\qquad\kappa(\textbf{A}\mathbf{M}^{-1})<<\kappa(\textbf{A}).$

Besides the computation of the matrix $\mathbf{M}^{-1}$, preconditioners
generally improve the numerical stability of direct and iterative solvers; for
iterative solvers, preconditioners generally increase the number of operations
performed for each iteration but reduce the number of iterations necessary to
converge.

###  2.2 Case studies

As case study, we consider the Laplace equation $\Delta u=0$ on a 2D/3D domain
$\Omega$ with Dirichlet condition $u=f$ on the boundary of $\Omega$. If the
input discrete domain is represented as a regular grid, then the Laplace
equation is discretised with a finite difference scheme and the corresponding
coefficient matrix is sparse and has a regular structure. If the input domain
is discretised as an irregular grid (e.g., a tetrahedral mesh as in Fig.
2(a)), then we apply a finite element discretisation [33]; the corresponding
coefficient matrix is still sparse but it has an irregular structure. For
tetrahedral meshes, the Laplacian matrix
$\mathbf{A}:=\mathbf{B}^{-1}\mathbf{L}$ is defined as the product between the
diagonal matrix $\mathbf{B}$, whose entries are the sum of the volume of the
tetrahedra incident at a given vertex, and $\mathbf{L}$ is the Laplacian
matrix with entries
$L(i,j):=w(i,j):=\frac{1}{6}\sum_{k=1}^{n}l_{k}\cot\alpha_{k}$ for each edge
$(i,j)$, $L(i,i):=-\sum_{j\in N(i)}w(i,j)$, and zero otherwise [31]. For
regular grids, $w(i,j):=1$. In all the aforementioned cases, the Laplace-
Beltrami operator is discretised as the $n\times n$ _Laplacian matrix_

$A(i,j):=\left\\{\begin{array}[]{ll}w(i,j)&(i,j)\textrm{ edge};\\\
-\sum_{k}w(i,k)&i=j,\end{array}\right.$

which is sparse and positive semi-definite, but not necessarily symmetric.
Then, the Laplace equation is discretised as the sparse linear system
$\textbf{A}\textbf{u}=\textbf{b}$, where u is the array of the values of the
solution at the nodes of the input grid and $\mathbf{b}$ is the right-hand
side term that defines the boundary condition. Finally, we do not make
assumptions on the coefficient matrix and/or constraints on the connectivity
of the input grid, such as a maximum degree of the nodes.

(a) | (b)
---|---
Figure 2: (a) Tetrahedralisation of a 3D domain and (b) sparsity pattern of
the corresponding Laplacian matrix.

Recalling that the _non-zero density_ of a matrix is defined as the percentage
of its non-zero entries with respect to the total number of elements of the
matrix and the _sparsity pattern_ as the structure of the non-zero entries in
the coefficient matrix, the Laplacian matrix on irregular grids (Table 1) has
an arbitrary sparsity pattern and a non-zero density percentage that is
generally higher than on regular grids, as a matter of the generally higher
number of edges incident to a given node. Indeed, the sparsity pattern of the
Laplacian matrix associated with irregular grids has an impact on the data
distribution among processes, MPI communications, and processes work load.
However, the reduced non-zero density percentage can be exploited to improve
the arithmetic density of the operations through an appropriate column
reordering. These considerations further motivate a comparison of solvers of
sparse linear systems associated with PDEs discretised on a regular/irregular
grid.

Table 1: Grid size and sparsity of the coefficient matrix of the input sparse
linear systems.

Regular grids
---
Domain |  $\sharp$Nodes |  $\sharp$ Mat. Rows |  $\sharp$ Non-zeros |  Non-zeros %
Cube 1 | $128\times 128\times 128$ | 2 097 152 | 14 099 408 | $3\text{\times}{10}^{-4}$
Cube 2 | $256\times 256\times 256$ | 16 777 216 | 115 099 600 | $4.1\text{\times}{10}^{-5}$
Cube 3 | $512\times 512\times 512$ | 134 217 728 | 930 123 728 | $5\text{\times}{10}^{-6}$
Irregular grids
Domain |  $\sharp$Nodes |  $\sharp$ Mat. Rows |  $\sharp$ Non-zeros |  Non-zeros %
Sphere 1 | 2 094 977 | 2 094 977 | 33 225 967 | $7\text{\times}{10}^{-4}$
Sphere 2 | 2 094 834 | 2 094 834 | 61 055 286 | $1.4\text{\times}{10}^{-3}$

###  2.3 Metrics and testbeds

To evaluate the solvers’ performances, we measure the following metrics

  * •

_execution time_ expressed in seconds;

  * •

_FLOPS_ \- FLoating point Operations Per Second - expressed in Giga;

  * •

_efficiency $(n)$_=$\frac{speedup(n)}{n}$ with $n$ processes;

  * •

_FLOPS efficiency $(n)$_$=\frac{F(n)}{F(1)}/n$ where $F(1)$ and $F(n)$ are
FLOPS with $1$ and $n$ processes, respectively;

In order to discuss the accuracy of direct and iterative solvers, we measure
the approximation error between the exact and the computed solution as

$x_{err}=\frac{\left\|\mathbf{x}_{ground\\_truth}-\mathbf{x}_{computed}\right\|_{2}}{\left\|\mathbf{x}_{ground\\_truth}\right\|_{2}},$
(1)

where the discrete ground-truth solution of the input PDE is given by sampling
the analytic solution at the vertices of the input domain. As exit condition
of iterative solvers, we consider the relative error

$\frac{\left\|\mathbf{A}\mathbf{x}-\mathbf{b}\right\|_{2}}{\left\|\mathbf{b}\right\|_{2}}\leq\epsilon,$
(2)

where $\epsilon$ is a given threshold.

The memory requirement of each method is measured by the total memory (i.e.,
including support matrices and vectors) allocated in the RAM. The allocated
memory is expressed in MB. These metrics will be used for an analysis of the
performances (Sect. 3) and for a comparison (Sect. 4) of direct and iterative
solvers, in order to identify which is more suitable on the basis of the
application and the resources available. All performance tests have been run
on the CINECA cluster Marconi, based on Intel Xeon product family and
classified at the 12-th position on Top500 ranking [37] in November 2016. In
particular, we have considered the Broadwell partition, composed by 2x18-cores
Intel Xeon E5-2697 v4 (Broadwell) at 2.30 GHz, 1512 nodes, 36 cores/node.
Total cores: 54432; 128 GB/node of RAM.

Figure 3: SuperLU execution time profiled according to its five main
operations.

##  3 Experimental results

We discuss our experiments on direct (Sect. 3.1) and iterative (Sect. 3.2)
solvers.

###  3.1 Experimenting direct solvers

According to our preliminary comparison [11] between MUMPS [5, 6] and
SuperLU_DIST111Version SuperLU_DIST 5.2.0 is the parallel extension to the
serial SuperLU for parallel distributed memory architectures [27]. The latter
obtained better performance and has been selected for a further investigation
in this work, and for the sake of simplicity and readability since now on we
will refer it as SuperLU. It subdivides the computation in five main steps

  1. 1.

_data distribution_ among processes, which implements a cyclic block pattern
of different size, according to the non-zero patterns of the input coefficient
matrix;

  2. 2.

_column permutation_ , which reorders the columns of the coefficient matrix to
optimise the number of non-zero entries after factorisation;

  3. 3.

_symbolic factorisation_ , which provides the number and position of non-zero
entries of the ($\mathbf{L},\mathbf{U}$) matrices;

  4. 4.

_factorisation_ , which factorises the coefficient matrix in the LU form,
where $\mathbf{L}$ and $\mathbf{U}$ are a lower and upper triangular matrix,
respectively;

  5. 5.

_solver_ , which computes the solution of the triangular systems.

As input data sets (Table 1), we consider a regularly tessellated cube and an
irregularly tessellated sphere. According to Table 2, the LU factorisation
increases the fill-in ratio of the input coefficient matrix by 2 and 3 orders
of magnitude for regular and irregular domains, respectively. This phenomenon
leads to higher memory and computational costs for the solution of the
corresponding linear systems.

Table 2: Size, non-zeros density and fill-in percentage of the coefficient matrix and its LU factorisation. Domain | A rows |  A Non-zeros % | LU fill-in %
---|---|---|---
Cube 1 | 2097152 |  $3.2\text{\times}{10}^{-6}$ |  $9.6\text{\times}{10}^{-4}$
Sphere 1 | 2094977 |  $7.5\text{\times}{10}^{-6}$ |  $1.7\text{\times}{10}^{-3}$

In Fig. 3, we report the execution time achieved by varying the number of
processes up to 576 processes and profiled on the base of the aforementioned
operations. Due to memory requirements for the allocation of the factorised
matrices, it was not possible to solve the system with a single processor;
indeed, for the performance evaluation we compare the execution time on $p$
processors with respect to the execution time on $p_{1}$ processors,
$p_{1}>1$. Memory allocation for both direct and iterative solvers are
discussed in Sect. 4.1.

The system resolution on a regular grid requires more than 20 minutes with 2
processes, reduces to less than 1 and a half minute with 72 processes, and 45
seconds with 576 processes. For irregular grids, the numerical solver requires
more than 50 minutes with 2 processes, reduces to 2 minutes and a half with 72
processes, and 37 seconds with 576 processes. Growing from 288 to 576
processes has a little impact on time reduction, especially for regular grids.
This aspect is partially due to the sequential executions of the _column
permutation_ and of the _symbolic factorisation_. Since these operations do
not scale on the processes, their impact on the global computation becomes
heavier while increasing the number of processes. On regular grids, the
aggregated impact of both operations is 3% with 2 processes, and it increases
to 55% with 576 processes; on irregular grids, it grows from 1.7% with 2
processes to 62% with 576 processes. On the contrary, _data distribution_ and
_solver_ smoothly scale and have a negligible effect. In fact, the aggregated
impact counts less than 3% on the global time. Indeed, the _factorisation_
dominates algorithm’ complexity, [27, 26, 29].

To further investigate the _factorisation_ operation (yellow parts of Fig. 3),
in Table 3 we report FLOPS and efficiency on regular and irregular grids
obtained by varying the number of processes. We note that the metric is
normalised with 2 processes, as memory allocation problems occurred with 1
process. The values increase from about 32 GFLOPS on 2 processes to about 1
TFLOPS (11.1% efficiency) and to 2.7 TFLOPS on 576 processes (30% efficiency)
on regular and irregular domains, respectively. It follows that the FLOPS
metrics is interesting even if the efficiency is limited, since reaching the
Tera scale is a plus. In general, better results are achieved on irregularly
tessellated domains, as they have an improved arithmetic density for direct
solvers, while maintaining limited the impact of the memory pattern and
communication overhead.

Table 3: Direct solvers, FLOPS, and factorisation efficiency. Processes | Regular grid | Irregular grid
---|---|---
| GFLOPS | Efficiency [%] | GFLOPS | Efficiency [%]
1 | - | - | - | -
2 | 32 | 100 | 32 | 100
4 | 64 | 99.8 | 66 | 100
8 | 119 | 93.2 | 129 | 99.6
18 | 189 | 65.6 | 272 | 92.8
36 | 288 | 49.8 | 504 | 85.9
72 | 467 | 40.4 | 758 | 64.6
144 | 670 | 29.0 | 1330 | 56.6
288 | 903 | 19.5 | 2073 | 44.2
576 | 1033 | 11.1 | 2765 | 29.4

###  3.2 Experimenting iterative solvers

We selected PETSc [9] to test 5 iterative solvers (BICGSTAB, GMRES, IBICGSTAB,
TFQMR, CGR) of sparse linear systems without applying preconditioners and
combined with 3 preconditioners (Hypre, ASM, Block-Jacobi) [34]. According to
the results in [11], the combination of the iterative BICGSTAB method, which
is based on Lanczos bi-orthogonalisation, with the Block-Jacobi preconditioner
has been selected as solver to be further investigated, as it provided better
performance figures. We briefly recall that Block-Jacobi is a general purpose
preconditioner, which reduces the condition number of the coefficient matrix
without optimising its sparsity pattern and without affecting the non-zero
density. Indeed, we test the solver on the _Cube 1_ and _Spere 1_ domains, and
consider its convergence to a solution with $\epsilon=10^{-12}$ in Eq. (2). In
Fig. 4, we report the execution time of the BICGSTAB algorithm profiled on the
basis of the main matrix and vector operations. Execution time is very low on
both regular and irregular grids, however, the algorithm has a different
behaviour on regular and irregular grids. In the first case, the algorithm
requires about 17 seconds to converge to a solution with 1 process, and less
than half second with 576 processes; this smooth scalability can be improved
on higher matrix dimension, i.e. larger data sets, as discussed in Sect. 4.1.
On irregular grids, BICGSTAB requires about 2 and half minutes to converge
with 1 process, reduced to 13 seconds with 144 process; however, with 576
increases up to about 3 minutes. In fact, some of the main operations (e.g.,
_Matmult_ and _Vecdot_ , a matrix-vector multiplication and a vector scalar
product, blue and light-blue parts in Fig. 4) do not properly scale when
increasing the number of processes.

Figure 4: BICGSTAB execution time profiled with matrix and vector operations.

In Table 4, we report the performances of BICGSTAB on regular and irregular
grids in terms of FLOPS, related efficiency, and number of iterations; on
regular grids, the algorithm has an efficiency of 35% (13.5 GFLOPS) with 36
processes and of 17% (106 GFLOPS) with 576 processes. On irregular grids, the
trend is significantly worst, with an efficiency of 17% (2.5 GFLOPS) with 36
processes and lower than 1% (and 0.5 GFLOPS) with 576 processes. Furthermore,
the number of iterations needed to the BICGSTAB algorithm to converge grows as
we increase the number of processes; this aspect has an impact on the solving
time and, consequently, on the algorithm efficiency.

Table 4: FLOPS, efficiency, and number of iterations of BICGSTAB on regular
and irregular grids.

|  Regular grid \- Cube 1 |  Irregular grid \- Sphere 1
---|---|---
Processes | GFLOPS | Iterations | FLOPS Eff. [%] | GFLOPS | Iterations | FLOPS Eff.
1 | 1.1 | 128 | 100 | 0.4 | 195 | 100
2 | 2.4 | 152 | 108 | 0.8 | 275 | 99
4 | 4.3 | 169 | 98 | 1.3 | 463 | 80
8 | 6.3 | 168 | 72 | 1.6 | 504 | 50
18 | 6.8 | 164 | 35 | 1.7 | 546 | 23
36 | 13.2 | 183 | 34 | 2.5 | 574 | 17
72 | 21.6 | 205 | 28 | 4.2 | 528 | 14
144 | 45.2 | 242 | 29 | 8.0 | 591 | 14
288 | 92.2 | 231 | 30 | 3.1 | 521 | 3
576 | 104.0 | 270 | 17 | 0.5 | 539 | 1

(a) Regular grids

(b) Irregular grids

Figure 5: Matmult communications and messages’ length for regular/irregular
grids.

Since the lower regularity of the sparsity pattern of the coefficient matrix
has a deep impact on the slow convergence of the solver, we analyse the trend
of the main operations. In Fig. 5, we propose an example of the communication
figures in term of MPI message number and average message length for the
_Matmult_ operation. Increasing the number of processes, the number of
communications increases while their size reduces. However, on regular grids
(a) the gap is closer than on irregular grids (b), In fact, the former message
length varies from an average of $100$ KB to $10$ KB and the message number
varies from $10^{3}$ to $10^{6}$, as for the latter message length varies from
an average of $1$ MB to $1$ KB and message number varies from $10^{3}$ to
$10^{8}$. The overhead paid to send a high number of short messages (Fig.
5(b), last column) is not balanced by a more intensive computational
counterpart neither by a more accurate solution leading to a faster
convergence. Actually, this trend is confirmed also on more demanding linear
systems, i.e. on _Sphere 2_.

##  4 Computational performances and numerical accuracy

Besides scalability and performance figures, we propose additional experiments
to characterise the properties of linear and iterative solvers peculiarities
in terms of memory requirements (Sect. 4.1), approximation error (Sect. 4.2),
and multiple right-hand side terms (Sect. 4.3). We also discuss how the choice
of a polynomial instead of a linear basis of the FEM method affects the
performances of iterative solvers.

###  4.1 Memory requirements: regular versus irregular grids

We analyse memory requirements of linear solvers, as they represent a
bottleneck for the execution for direct solvers and for the scalability of
iterative solvers.

Table 5: Memory requirements (expressed in MB) of direct and iterative solvers on the regular grid _Cube 1_ and the irregular grid _Sphere 1_ \- average per process.  Processes | Regular grid | Irregular grid
---|---|---
| Direct [MB] | Iterative [MB] | Direct [MB] | Iterative[MB]
1 | - | 588 | - | 1020
2 | 17976 | 340 | 37603 | 529
4 | 9284 | 196 | 17512 | 298
8 | 4899 | 123 | 8471 | 192
18 | 2405 | 82 | 3972 | 134
36 | 1393 | 66 | 1850 | 112
72 | 831 | 58 | 1198 | 76
144 | 580 | 53 | 694 | 61
288 | 417 | 50 | 428 | 54
576 | 369 | 49 | 314 | 53

Table 5 reports the allocated memory of direct and iterative solvers for _Cube
1_ and _Sphere 1_ , expressed as the average memory allocated per process.
Both methods have a good distribution ratio of data among processes; average
and maximum values are quite close, apart from few cases for direct solvers.
However, direct solvers allocate one or two orders of magnitude more than
iterative solvers, as a very large amount of memory is allocated for the
factorised matrices; for instance, direct solvers with two processes need 17.5
GB of memory and 37 GB versus 340 MB and 530 MB required by iterative solvers
on regular and irregular grids, respectively. Indeed, in the latter case both
methods double the memory costs, without pushing to a higher order.

Table 6: Scalability results of iterative solvers when increasing the
dimension of the coefficient matrix, with Time, GigaFLOPS (GF.) and FLOPS
efficiency (F. Eff.) metrics.

Domain | Cube 1 | Cube 2 | Cube 3
---|---|---|---
Procs | Time[s] | GF. | F. Eff.[%] | Time[s] | GF. | F. Eff.[%] | Time[s] | GF. | F. Eff.[%]
1 | 17.4 | 1 | 100 | 1047.1 | 0.3 | 100 | 4361.8 | 0.3 | 100
144 | 0.73 | 46 | 29 | 12.2 | 34 | 72 | 141 | 41 | 89
288 | 0.32 | 94 | 30 | 9.7 | 55 | 57 | 123 | 58 | 58
576 | 0.34 | 106 | 17 | 6.7 | 81 | 42 | 73.0 | 110 | 54

From these tests, it follows that iterative solvers definitely overwhelm
direct solvers in terms of memory cost, and represent the unique possible
choice to manage large coefficient matrices. In fact (Sect. 3.1), direct
solvers have incurred in memory limits for the computation of the LU
factorisation when largely increasing the number of nodes of the input grid.
On the contrary, iterative solvers generally do not suffer of memory limits,
once the coefficient matrix has been allocated.

We further test iterative solvers considering larger regular domains.
According to the results in Table 6, the scalability of the algorithm improves
with the dimension of the coefficient matrix: BICGSTAB requires (i) 18 minutes
with 1 process and 7 seconds with 576 processes on _Cube 2_ with a $256^{3}$
coefficient matrix; (ii) about 1 hour and 13 minutes with 1 process and about
1 minute on _Cube 3_ with a $512^{3}$ coefficient matrix. For FLOPS, tests
outlined 17% efficiency with _Cube 1_ , 42% with _Cube 2_ and 54% with _Cube
3_ when considering 576 processes. This result derives from an improved
balancing between computation and communications: Fig. 6 provides the
percentage of time spent in MPI communications and computation, while varying
the number of processes and the size of the coefficient matrix. Better ratios
are achieved with more demanding data sets and a variable number of processes,
as reflected in FLOPS efficiency values.

Figure 6: Timings of the computation and communication phases by varying the
number of processes and the size of the coefficient matrix; tests consider
regular grids.

###  4.2 Error analysis

To study how the number of iterations and the stopping criteria influence the
accuracy and the execution time of iterative solvers, we set the relative
error for the stopping condition in Eq. (2) equal to
$\epsilon:=10^{-8},10^{-12},10^{-15}$. For direct solvers based on the
SuperLU, we consider the error metrics (1) with $x_{err}\approx 10^{-15}$. Our
tests are performed with 72 parallel processes, i.e., 2 node of Marconi
Supercomputer and the _Cube 1_ domain.

Fig. 7 compares the approximation accuracy of different iterative solvers
without preconditioning; in particular, we consider BICGSTAB, GMRES, and IDR
solvers. This comparison is intended to understand if the selection of a
specific iterative algorithm may have an impact on the accuracy. Here, the
horizontal line represents the approximation error $x_{err}$ of SuperLU.
Iterative solvers may differ of one order in the approximation error and
BICGSTAB generally shows better results, which are further investigated by
considering additional preconditioners.

Figure 7: Solution accuracy ($y-$axis) of different iterative solvers without
preconditioning, with respect to different values of $\epsilon$ ($x$-axis).

Fig. 8 compares the approximation accuracy of BICGSTAB with the Block-Jacobi
and Hypre preconditioners and without preconditioners, by varying $\epsilon$.
Hypre needs a larger time (about 25 seconds) to reach a very accurate
solution, due to the strong preconditioning phase, while Block-Jacobi has a
good approximation accuracy against a negligible execution time (0.01
seconds), which leads to a $50\%$ reduction of the number of iterations with
respect to the solution without preconditioning; for sake of readability in
the Figure we depict only the latter. The red bars show the solving time of
_Block Jacobi - BICGSTAB_ for different $\epsilon$; varying $\epsilon$ from
$10^{-8}$ to $10^{-15}$, the error reaches the order of direct solvers paid
with a little impact on execution time, which is still not comparable with
respect to direct solvers and provides advantages from memory point of view.
We further investigate Block Jacobi - BICGSTAB with respect to the
approximation accuracy.

Figure 8: ($y-$axis) Solution accuracy and timing of the BICGSTAB solver with
different preconditioners and values of $\epsilon$ ($x$-axis).

#### Further options: polynomial basis

An important parameter in the modelling of the applicative problem is the
degree of the polynomial basis used for the FEM discretisation of the input
PDE, which affects the sparsity pattern of the coefficient matrix
$\mathbf{A}$, and thus the solution of the sparse linear system. We analyse
the impact on the approximation accuracy when increasing the degree of the
polynomial basis from degree 1, i.e. the value used in all tests of this work,
to degree 4. Table 7 reports the impact for each degree, the coefficient
matrix $\mathbf{A}$ size and its non-zero numbers increase when passing from
polynomial degree 1 to degree 4, together with an increase of the number of
iterations and execution time ($\epsilon$ fixed to $10^{-12}$). The
approximation accuracy improves significantly from degree 1 to 2 slightly
affecting execution time, while approximation accuracy remains almost
unchanged with degree 3 and 4 with a notable impact on execution time.

Table 7: Block Jacobi - BICGSTAB results with different polynomial degrees, in
terms of execution time and accuracy ($\epsilon=10^{-12}$).

Basis | Nodes | Elem. | Non-zeros |  Solver [$s$] | Iter. | Sol. accur.
---|---|---|---|---|---|---
P1 | 86 425 | 86 425 | 602 969 | 0.11 | 358 |  $5.0\text{\times}{10}^{-6}$
P2 | 86 425 | 344 697 | 3 894 474 | 0.52 | 689 |  $2.3\text{\times}{10}^{-10}$
P3 | 86 425 | 774 817 | 13 148 592 | 3.32 | 938 |  $6.0\text{\times}{10}^{-10}$
P4 | 86 425 | 1 376 785 | 32 319 425 | 10.21 | 1260 |  $1.7\text{\times}{10}^{-9}$

The choice of the polynomial degree influences the approximation accuracy of
the solution; it could be selected on the basis of application requirements
but it is necessary to take into account the increased computational cost.

###  4.3 Handling multiple right-hand side terms

Different physical phenomenon can be described through PDEs discretised by
systems with the same coefficient matrix and multiple right-hand side (r.h.s.)
terms. For instance, we can consider the approximation of the equilibrium
configuration of an isotropic membrane, whose boundary changes its geometry
with respect to time. Therefore, we evaluate the execution time of direct and
iterative solvers of large sparse linear systems with multiple right-hand side
terms. Direct solvers become valuable in this case, since most of the
operations (e.g., LU factorisation and column permutation of the coefficient
matrix) are performed only once while the solution is computed for each r.h.s.
term with a backward and forward substitution in linear time. On the contrary,
iterative methods solve a new linear system for every right-hand side term,
without preserving any computational phase although the coefficient matrix
does not change. During the iterations, only the preconditioner is preserved;
according to Sect. 3.2, the negligible or predominant computation of the
preconditioner depends on the selection of a weak or a strong preconditioner.
SuperLU provides the following options

$\left\\{\begin{array}[]{ll}\mathbf{A}\mathbf{x}_{i}=\mathbf{b}_{i},\,i=1,\ldots,t&\textrm{(single
r.h.s terms)};\\\ \\\
\mathbf{A}\mathbf{X}=\mathbf{B},\quad\mathbf{X}:=[\mathbf{x}_{1},\ldots,\mathbf{x}_{t}],\quad\mathbf{B}:=[\mathbf{b}_{1},\dots,\mathbf{b}_{t}]&\textrm{(block
of r.h.s terms)}.\end{array}\right.$ (3)

The first equation corresponds to the case where r.h.s. terms is passed as one
vector per time; the second equation supposes that all the r.h.s. vectors are
known a-priori (e.g., they do not depend on the solution computed for the
previous r.h.s. terms) and that are passed as a _block_ of r.h.s. terms; i.e.,
a matrix $\mathbf{B}$. Generally, the matrix $\mathbf{B}$ is dense and non
symmetric; therefore, memory limits may be experienced.

For iterative solvers, we choose $\epsilon=10^{-8},10^{-12},10^{-15}$ on
regular grids and we choose only $\epsilon=10^{-8}$ on irregular grids. Tests
have been executed considering 2 nodes of Marconi, i.e. 72 parallel processes,
on the data-sets _Cube 1_ and _Sphere 1_. As performance metric, we propose
execution time varying the number of r.h.s. terms; this means that: (i) the
lines with a slower slope correspond to faster solvers, (ii) possible
intersections between two lines correspond to a change of behaviour, i.e. the
number of r.h.s. terms where one solver becomes more efficient (faster) than
the other.

We also estimate the number of terms such that direct and iterative solvers
require the same computation time. If $t$ is the number of right-hand side
terms, the computation time $T_{SLU}(t)$ of direct solvers of $t$ linear
systems is equal to the factorisation time $T_{SLU\\_FACT}$ of the coefficient
matrix $\mathbf{A}$ plus the solving time $T_{SLU\\_SOLVE}$ of the two
triangular systems for each r.h.s. term; i.e.,

$T_{SLU}(t)\approx(T_{SLU\\_FACT}+t\cdot T_{SLU\\_SOLVE}).$

Similarly, the computation time $T_{IS}(t)$ of iterative solvers of $t$ linear
systems is equal to the preconditioning time $T_{IS\\_PC}$ plus the
coefficient matrix and the solving time $T_{IS\\_SOLVE}$ for each r.h.s. term;
i.e.,

$T_{IS}(t)\approx(T_{IS\\_PC}+t\cdot T_{IS\\_SOLVE}).$

Then, the number $t_{0}$ of r.h.s. terms such that the execution time of
linear and iterative solvers is comparable is equal to

$t_{0}\approx\frac{T_{SLU\\_FACT}-T_{IS\\_PC}}{T_{IS\\_SOLVE}-T_{SLU\\_SOLVE}}.$

Table 8 reports the value $t_{0}$ for tests performed on regular/irregular
grids and by varying $\epsilon$. Figs. 9, 10 depict results on
regular/irregular grids from direct (two options) and iterative (two
preconditioned) solvers. As for regular grids, direct solvers applied to
linear systems with multiple r.h.s. terms outperform direct solvers applied to
a r.h.s vector per time. BICGSTAB preconditioned with Hypre has the worst
performance, since the coefficient matrix is already regular and well
distributed. If the exit condition (i.e., the parameter $\epsilon$) is
relaxed, then iterative solvers have a good performance on regular grids and
these results are better or comparable with respect to direct solvers; when
$\epsilon$ is reduced, the situation changes. In fact, with
$\epsilon=10^{-8}$, BICGSTAB preconditioned with Block Jacobi (blue line) has
a smaller slope than direct method with single r.h.s. terms (red line); direct
solvers with a block of r.h.s. terms (violet line) are faster considering more
than 290 r.h.s. terms. BICGSTAB preconditioned with Hypre (yellow line) has a
larger execution time; indeed, direct solvers obtain more effective results
after 210 and 100 r.h.s. terms for the single and block r.h.s. terms,
respectively.

Table 8: Average number of r.h.s. terms such that iterative and direct solvers
have a comparable execution time; Inf and Sup indicate that the method on the
horizontal line has always the best/worst performance with respect to vertical
intersection.

| Regular grid | Irregular grid
---|---|---
| $\epsilon=1e-8$ | $\epsilon=1e-12$ | $\epsilon=1e-15$ | $\epsilon=1e-8$
|  BJacobi |  Hypre |  SLU-single |  SLU-multiple |  BJacobi |  Hypre |  SLU-single |  SLU-multiple |  BJacobi |  Hypre |  SLU-single |  SLU-multiple |  BJacobi |  Hypre |  SLU-single |  SLU-multiple
BJacobi |  - |  inf |  inf |  290 |  - |  inf |  680 |  165 |  - |  inf |  260 |  120 |  - |  41 |  16 |  15
Hypre |  sup |  - |  210 |  100 |  sup |  - |  175 |  90 |  sup |  - |  125 |  75 |  41 |  - |  9 |  8
SLU-single |  sup |  210 |  - |  sup |  680 |  175 |  - |  sup |  260 |  125 |  - |  sup |  16 |  9 |  - |  sup
SLU-multiple |  290 |  100 |  inf |  - |  165 |  90 |  inf |  - |  120 |  75 |  inf |  - |  15 |  8 |  inf |  -

Increasing the solution accuracy ($\epsilon=10^{-15}$), iterative solvers
become less competitive; direct method with a block of r.h.s. terms provides
the best results with respect to the BICGSTAB method preconditioned with Block
Jacobi and with Hypre with more than 120 and 75 r.h.s. terms, respectively.
BICGSTAB preconditioned with Block Jacobi performs better than direct solvers
with single r.h.s. terms up to 260 vectors, while BICGSTAB preconditioned with
Hypre is interesting up to 125 r.h.s. terms. Intermediate $N$ values are
obtained for $\epsilon=10^{-12}$; results are reported in Table 8.

Figure 9: Comparison among several approaches in terms of execution time
($y-$axis) with a different number of r.h.s. terms ($x-$axis) on regular
grids.

Our tests (Fig. 10) on an irregularly tessellated domain and with
$\epsilon=10^{-8}$ lead to different results. Iterative solvers result less
competitive even with low $\epsilon$ values on irregular grids. BICGSTAB
preconditioned with Block Jacobi (blue line) and with Hypre (yellow line) are
effective only with a very low number r.h.s. terms; i.e. 16 and 8,
respectively. This result is due to their poor scalability figures and to the
very slow convergence on irregular grids, as already described in Sect. 3.2.
If we select an iterative solver, then a strong preconditioner (i.e., Hypre)
provides good results in terms of execution time when more than 40 r.h.s.
terms are selected. On the basis of these results, we do not further decrease
the value of $\epsilon$ on irregular domain.

##  5 Discussion

#### Direct versus iterative solvers

SuperLU has good efficiency results on both regular and irregular grids, with
an improvement on the latter. On regular grids, SuperLU has a shorter
computing timing due to a lower number of non-zero entries. On irregular
grids, SuperLU has better scalability results; in fact, a higher non-zero
density improves the arithmetic density. The irregular sparsity pattern is
managed by a column reordering, which allows a very good distribution of the
data. The actual bottleneck of SuperLU is represented by memory requirements
that are due to the LU factorisation; therefore, architectural limits can be
experienced. This situation may limit the exploitation of SuperLU in terms of
the size of the problem/application that can be solved, and in terms of memory
capacity of the underlying computational resources. However, SuperLU is
particularly efficient for the solution of large sparse linear systems with
multiple r.h.s. terms, and the approximation accuracy is close to the round-
off error, when considering well-conditioned coefficient matrices.

Figure 10: Comparison among several approaches in terms of time ($y-$axis)
with a different number of r.h.s. terms ($x-$axis) on irregular grids.

As for iterative solvers, better performances are obtained on regular grids,
while irregular grids experience an increased number of iterations, due to a
more irregular sparsity pattern of the coefficient matrix. We observe that the
main operations of the solver do not scale when increasing the number of
processes; this aspect has an impact on the solving time and, consequently, on
the algorithm efficiency. Beside that, execution time of iterative solvers is
considerably lower with respect to direct method ones and memory limits never
occurred during our tests. The potential drawback of iterative solvers is a
generally lower numerical accuracy, which depends on the number of iterations
and the threshold selected for the stopping condition in Eq. (2). This
situation can be limited by the means of a strong preconditioner and/or a
higher polynomial basis degree in the discretisation of Laplace equation, both
paid with an increased computational cost and execution time.

Table 9: Main characteristics of direct and iterative linear solvers. |  Accuracy |  Execution time |  Memory requirements |  Data set |  Number r.h.s. terms | Architecture
---|---|---|---|---|---|---
|  |  |  |  |  One |  More |  Embedded |  Medium |  Large
Direct Solver |  High numerical accuracy
The approximation accuracy is close to the round-off error, in case of well-conditioned or preconditioned coefficient matrices. |  Expensive LU factorisation
Fast linear solver

Execution time is competitive with respect to other direct or iterative
solvers in case of multiple r.h.s. terms.

High scalability on regular and irregular grids |  High memory overload for the storage of the $(\mathbf{L},\mathbf{U})$ factors. |  Regular Grid

Irregular Grid |  |  ✓

✓ |  |  ✓

✓ |  ✓

✓
Iterative Solver |  Low numerical accuracy but tuneable according to stopping criteria (e.g., number of iterations, numerical accuracy).
|  Low numerical accuracy
The method is valuable when solving a single system.

In case of $t$ r.h.s. terms, the computational cost increases linearly with
respect to $t$.

Higher execution time in case of irregular grids. |  Low memory overload, once the coefficient matrix has been stored. Best choice in case of very large sparse linear systems, whose LU factorisation does not fit into the main memory. |  Regular Grid

Irregular Grid

|  ✓

✓ |  |  ✓

✓ |  ✓

✓ |  ✓

#### Moving towards multiple r.h.s. terms

In case of multiple r.h.s. terms, direct solvers generally outperform
iterative solvers even when the threshold on accuracy is relatively low.
Indeed, iterative solvers are more effective than direct solvers, in case of a
low number of r.s.h. terms. If the number of r.s.h. terms grows, direct
solvers are more efficient than iterative solvers. In Sect. 4.3, we have
identified the average number of terms that make the computation time of
direct and iterative solvers comparable; in this case, we have considered each
operation of both methods (i.e., factorisation, preconditioning) as a black
box, by measuring the execution time without making any assumption on the
performed operations.

We can further analyse this aspect from a computational point of view, i.e.,
by identifying the average number of r.h.s. terms such that direct solvers
have a computational cost comparable with respect to iterative solvers. In
case of $t$ r.h.s. terms and an input coefficient matrix with $m$ non-zero
entries, we have that direct solvers take $\mathcal{O}(m)$-time for the LU
factorisation of the input coefficient matrix, $\mathcal{O}(n)$-time for the
solution of a linear system with one r.h.s. term through the backward and
forward substitution. Indeed, the overall computational cost is
$\mathcal{C}_{ds}:=\mathcal{O}(m+tn)$. Iterative solvers solve one linear
systems for each r.h.s. term; indeed, the overall computational cost is
$\mathcal{O}(\sum_{i=1}^{t}k(i)n)$, where $k(i)$ is the number of iterations
for the solution of the linear system with the $i$-th r.h.s. term. Indicating
with $k_{max}$ the maximum number of iterations, the computational cost is
$\mathcal{C}_{is}:=\mathcal{O}(k_{\max}tn)$.

The computational cost of direct solvers is lower than iterative solvers if
$\mathcal{C}_{is}>\mathcal{C}_{ds}$, i.e., $t>t_{0}:=m/(n(k_{\max}-1))$.
Otherwise, direct solvers are computationally more expensive than iterative
solvers. If the iterative solver is slowly converging (i.e., $k_{\max}$ is
large) and the LU factorisation fits into the main memory, then direct solvers
are generally more efficient. This remark and the results of our experiments
(Sect. 3) confirm the importance of the preconditioning of the coefficient
matrix for iterative solvers, which will affect the value $t_{0}$ of a factor
$\varphi(n)/(k_{max}-1)$, where $\varphi(n)$ is the cost of the
preconditioning step (e.g., $\varphi(n)=\mathcal{O}(n)$, for most of the
preconditioners).

#### Tools usability

Another interesting aspect is the usability of the tools in terms of
developer’ skills and programming effort. SuperLU and PETSc manage the
parallelism in a transparent manner with respect to the users, and provide
statistics for a high-level profiling of the performances. Since both tools
offer the possibility to solve sparse linear systems with few lines of code,
their easiness of use is comparable. PETSc represents a more general framework
to solve large sparse linear systems and it can be integrated with different
external tools, including SuperLU. These considerations imply a larger
community of users, a more detailed documentation available online, and
support forums. Finally, the performances of SuperLU, used as a stand alone
tool or through PETSc interface, are similar.

#### Rule of thumbs

We summarise these considerations in Table 9 (pag. 9), which presents a
general rule of thumbs aimed at supporting application specialists. It
provides an empirical base to understand the main peculiarities of direct and
iterative solvers, thus driving the effort of application specialists towards
the selection of the solver that better fits their requirements.

Taking into account the application requirements (e.g., number of r.h.s.
terms, coefficient matrix characteristics, target approximation accuracy)
together with the underlying computational resources (e.g., memory capacity),
application specialists have a first insight on the possible options. If the
PDEs have to be run on embedded systems (i.e., a resource with limited
capabilities), then iterative solvers should be considered as first option. If
the input problem requires a high accuracy, then direct solvers should be
preferred. If the problem leads to a very large coefficient matrix, then
iterative solvers should be considered. In case of multiple r.h.s. terms,
previous considerations and in Sect. 4.3 guide the user in the selection of
the solver with the lowest computation time.

##  6 Conclusions

The solution of PDEs through approximate techniques typically results in the
solution of large sparse linear systems. In this context, we selected the
Laplace equation on a 2D/3D domain as a standard problem for several physical
applications, and we proposed a study about widespread parallel tools for the
solution of sparse linear systems with direct and iterative solver. We
evaluated their performances by taking into account specific peculiarities of
the application domain, e.g. the discretisation of the problem on regular and
irregular input domain, the sparsity of the coefficient matrices, etc. We have
further compared direct and iterative solvers in terms of memory requirements,
numerical accuracy, and their effectiveness when handling of multiple r.h.s.
terms.

The final aim of the work is to support application-domain specialists by
providing a first insight on how solvers respond to their applications without
the need of performing a large number of experimentation and tests. The choice
between direct/iterative solvers has to take into account the knowledge of the
application specialists with respect to the underlying problem, and can be
based on different elements, such as the properties of the coefficient matrix
(e.g., sparsity, symmetry, positive definiteness), single/multiple r.h.s.
terms, and matrix size. According to these considerations, we proposed a
general rule of thumbs targeted on production-oriented programmers, i.e.
application-domain specialists exploiting off-the-shelf tools provided by
optimisation-oriented programmers. Our methodological approach can be applied
with different classes of PDEs or with approximation schemes. To the best of
our knowledge, a systematic analysis and the definition of a rule of thumbs
have not been proposed by scientific community while provides a pragmatic and
convenient support in the everyday work.

As future work, we plan to explore the exploitation of GPUs-based tools for
the solution of sparse linear systems; we are also interested in improve the
understanding of the influence of the discretisation of the input PDE with
piecewise linear and polynomial basis functions. In this paper, we mainly
relates it with the concept of accuracy, while we would study the wider impact
on computational costs and performance.

## Acknowledgments

This work has been partially supported by the Programme “ _Advanced Studies in
Scientific Computing_ ”, University La Sapienza - Rome, and the H2020 ERC
Advanced Grant CHANGE, Contract 694515. Tests have been supported by CINECA
through the ISCRA-C Project HPC-PDE.

## References

  * [1] Abdelfattah, A., Ltaief, H., Keyes, D., and Dongarra, J. Performance optimization of sparse matrix-vector multiplication for multi-component pde-based applications using gpus.  Concurrency and Computation: Practice and Experience 28, 12 (2016), 3447–3465.
  * [2] Agullo, E., Buttari, A., Guermouche, A., and Lopez, F. Multifrontal qr factorization for multicore architectures over runtime systems.  In European Conference on Parallel Processing (2013), Springer, pp. 521–532.
  * [3] Agullo, E., Buttari, A., Guermouche, A., and Lopez, F. Implementing multifrontal sparse solvers for multicore architectures with sequential task flow runtime systems.  ACM Transactions on Mathematical Software (TOMS) 43, 2 (2016), 13\.
  * [4] Amestoy, P. R., Davis, T. A., and Duff, I. S. An approximate minimum degree ordering algorithm.  SIAM Journal on Matrix Analysis and Applications 17, 4 (1996), 886–905.
  * [5] Amestoy, P. R., Duff, I. S., Koster, J., and L’Excellent, J.-Y. A fully asynchronous multifrontal solver using distributed dynamic scheduling.  SIAM Journal on Matrix Analysis and Applications 23, 1 (2001), 15–41.
  * [6] Amestoy, P. R., Guermouche, A., L’Excellent, J.-Y., and Pralet, S. Hybrid scheduling for the parallel solution of linear systems.  Parallel Computing 32, 2 (2006), 136–156.
  * [7] Asanovic, K., Bodik, R., Demmel, J., Keaveny, T., Keutzer, K., Kubiatowicz, J., Morgan, N., Patterson, D., Sen, K., Wawrzynek, J., et al. A view of the parallel computing landscape.  Communications of the ACM 52, 10 (2009), 56–67.
  * [8] Ashcraft, C., Grimes, R., Pierce, D., Wah, D., and Works, B. P. Solving linear systems using spooles 2.2.  Boeing Phantom Works (2002).
  * [9] Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Kaushik, D., Knepley, M. G., May, D. A., McInnes, L. C., Mills, R. T., Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H. PETSc users manual.  http://www.mcs.anl.gov/petsc, 2018.
  * [10] Buttari, A., Dongarra, J., Kurzak, J., Langou, J., Luszczek, P., and Tomov, S. The impact of multicore on math software.  In International Workshop on Applied Parallel Computing (2006), Springer, pp. 1–10.
  * [11] Cammarasana, S., Carlini, E., Clematis, A., Galizia, A., and Patanè, G. High performance computing for the efficient solution of pdes on arbitrary domains.  Master’s thesis, Master in Scientific Calculus, Universitá La Sapienza di Roma, https://www.mat.uniroma1.it/sites/default/import-files/didattica/master/TESI/1617/Cammarasana.pdf, 2016-217.
  * [12] Carracciuolo, L., Casaburi, D., D’Amore, L., D’Avino, G., Maffettone, P., and Murli, A. Computational simulations of 3d large-scale time-dependent viscoelastic flows in high performance computing environment.  Journal of Non-Newtonian Fluid Mechanics 166, 23-24 (2011), 1382–1395.
  * [13] Chen, C., Pouransari, H., Rajamanickam, S., Boman, E. G., and Darve, E. A distributed-memory hierarchical solver for general sparse linear systems.  Parallel Computing 74 (2018), 49–64.
  * [14] Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., and Speleers, H. Robust and optimal multi-iterative techniques for iga collocation linear systems.  Computer Methods in Applied Mechanics and Engineering 284 (2015), 1120–1146.
  * [15] Dongarra, J., Abalenkovs, M., Abdelfattah, A., Gates, M., Haidar, A., Kurzak, J., Luszczek, P., Tomov, S., Yamazaki, I., and YarKhan, A. Parallel programming models for dense linear algebra on heterogeneous systems.  Supercomputing frontiers and innovations 2, 4 (2016), 67–86.
  * [16] Dongarra, J., Heroux, M. A., and Luszczek, P. High-performance conjugate-gradient benchmark: A new metric for ranking high-performance computing systems.  The International Journal of High Performance Computing Applications 30, 1 (2016), 3–10.
  * [17] Dongarra, J. J., and Walker, D. W. Software libraries for linear algebra computations on high performance computers.  SIAM review 37, 2 (1995), 151–180.
  * [18] Evans, L. C. Partial differential equations.  American Mathematical Society, Providence, R.I., 2010.
  * [19] Falgout, R. D., Jones, J. E., and Yang, U. M. The design and implementation of hypre, a library of parallel high performance preconditioners.  In Numerical solution of partial differential equations on parallel computers. Springer, 2006, pp. 267–294.
  * [20] Golub, G. H., and Van Loan, C. F. Matrix computations, vol. 3.  JHU Press, 2012.
  * [21] Hénon, P., Ramet, P., and Roman, J. Pastix: a high-performance parallel direct solver for sparse symmetric positive definite systems.  Parallel Computing 28, 2 (2002), 301–321.
  * [22] Heroux, M. A., Bartlett, R. A., Howle, V. E., Hoekstra, R. J., Hu, J. J., Kolda, T. G., Lehoucq, R. B., Long, K. R., Pawlowski, R. P., Phipps, E. T., et al. An overview of the trilinos project.  ACM Transactions on Mathematical Software (TOMS) 31, 3 (2005), 397–423.
  * [23] Karypis, G., and Kumar, V. A fast and high quality multilevel scheme for partitioning irregular graphs.  SIAM Journal on scientific Computing 20, 1 (1998), 359–392.
  * [24] Keller Tesser, R., Mello Schnorr, L., Legrand, A., Heinrich, F. C., Dupros, F., and Navaux, P. O. Performance modeling of a geophysics application to accelerate over-decomposition parameter tuning through simulation.  Concurrency and Computation: Practice and Experience (2017), e5012.
  * [25] Kreutzer, M., Thies, J., Röhrig-Zöllner, M., Pieper, A., Shahzad, F., Galgon, M., Basermann, A., Fehske, H., Hager, G., and Wellein, G. Ghost: building blocks for high performance sparse linear algebra on heterogeneous systems.  International Journal of Parallel Programming 45, 5 (2017), 1046–1072.
  * [26] Li, X. S. An overview of superlu: Algorithms, implementation, and user interface.  ACM Transactions on Mathematical Software (TOMS) 31, 3 (2005), 302–325.
  * [27] Li, X. S., and Demmel, J. W. SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems.  ACM Trans. Mathematical Software 29, 2 (June 2003), 110–140.
  * [28] Liu, W., Li, A., Hogg, J. D., Duff, I. S., and Vinter, B. Fast synchronization-free algorithms for parallel sparse triangular solves with multiple right-hand sides.  Concurrency and Computation: Practice and Experience 29, 21 (2017), e4244.
  * [29] Liu, Y., Jacquelin, M., Ghysels, P., and Li, X. S. Highly scalable distributed-memory sparse triangular solution algorithms.  In 2018 Proceedings of the Seventh SIAM Workshop on Combinatorial Scientific Computing (2018), SIAM, pp. 87–96.
  * [30] Naumov, M. Incomplete-lu and cholesky preconditioned iterative methods using cusparse and cublas.  Nvidia white paper (2011).
  * [31] Patané, G. Star-laplacian spectral kernels and distances for geometry processing and shape analysis.  In Computer Graphics Forum (2016), vol. 35, Wiley Online Library, pp. 599–624.
  * [32] Puzyrev, V., Koric, S., and Wilkin, S. Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems.  Computers & Geosciences 89 (2016), 79–87.
  * [33] Quarteroni, A. Numerical models for differential problems, vol. 16.  Springer, 2017.
  * [34] Saad, Y. Iterative methods for sparse linear systems, vol. 82.  siam, 2003.
  * [35] Schenk, O., and Gärtner, K. Solving unsymmetric sparse systems of linear equations with pardiso.  Future Generation Computer Systems 20, 3 (2004), 475–487.
  * [36] Song, F., and Dongarra, J. A scalable approach to solving dense linear algebra problems on hybrid cpu-gpu systems.  Concurrency and Computation: Practice and Experience 27, 14 (2015), 3702–3723.
  * [37] Strohmaier, E. Top500 supercomputer.  In Proceedings of the 2006 ACM/IEEE conference on Supercomputing (2006), ACM, p. 18.
  * [38] Turek, S., Göddeke, D., Becker, C., Buijssen, S. H., and Wobker, H. Feast—realization of hardware-oriented numerics for hpc simulations with finite elements.  Concurrency and Computation: Practice and Experience 22, 16 (2010), 2247–2265.
  * [39] Wang, S., Li, X. S., Rouet, F.-H., Xia, J., and De Hoop, M. V. A parallel geometric multifrontal solver using hierarchically semiseparable structure.  ACM Transactions on Mathematical Software (TOMS) 42, 3 (2016), 21\.


# AIRCADE: an Anechoic and IR Convolution-based Auralization Data-compilation
Ensemble

Túlio Chiodi
School of Acoustical Engineering
Federal University of Santa Maria, Brazil
<EMAIL_ADDRESS>
& Arthur Nicholas dos Santos
School of Electrical and Computer Engineering
University of Campinas, Brazil
<EMAIL_ADDRESS>
& Pedro Martins
School of Electrical and Computer Engineering
University of Campinas, Brazil
<EMAIL_ADDRESS>
& Bruno Sanches Masiero
School of Electrical and Computer Engineering
University of Campinas, Brazil
<EMAIL_ADDRESS>

###### Abstract

In this paper111This work was partially supported by the São Paulo Research
Foundation (FAPESP), grants $\\#$2017/08120-6 and $\\#$2019/22795-1., we
introduce a data-compilation ensemble, primarily intended to serve as a
resource for researchers in the field of dereverberation, particularly for
data-driven approaches. It comprises speech and song samples, together with
acoustic guitar sounds, with original annotations pertinent to emotion
recognition and Music Information Retrieval (MIR). Moreover, it includes a
selection of impulse response (IR) samples with varying Reverberation Time
(RT) values, providing a wide range of conditions for evaluation. This data-
compilation can be used together with provided Python scripts, for generating
auralized data ensembles in different sizes: tiny, small, medium and large.
Additionally, the provided metadata annotations also allow for further
analysis and investigation of the performance of dereverberation algorithms
under different conditions. All data is licensed under Creative Commons
Attribution 4.0 International License.

_Keywords_ Speech Emotion Recognition $\cdot$ Song Emotion Recognition $\cdot$
Music Information Retrieval $\cdot$ Auralization $\cdot$ Dereverberation

##  1 Introduction

Reverberation is the persistence of sound in a space after the source ceases
to emit it. It is mainly caused by reflections of the sound waves off the
surfaces within the space and gradually dissipates over time. Since the
characteristics of the space in which reverberation occurs can influence the
duration and intensity of this phenomenon, the Reverberation Time (RT) is
considered an important aspect of room acoustics (Beranek (2004); Lapointe
(2010)).

In the context of entertainment, reverberation can affect the Human Auditory
System (HAS) in positive ways, since some degree of it can help to enhance the
perceived loudness and richness of sounds. This occurs because the multiple
reflections can create a sensation of spaciousness and immersion, which can be
aesthetically pleasing, particularly in music. On the other hand, in the
context of communications, excessive reverberation can affect the quality and
intelligibility of speech because the multiple reflections can create a
"smearing" effect that can make it difficult to distinguish individual sounds
and syllables. Moreover, the reduction of the overall Signal-to-Noise Ratio
(SNR) of a given sound source can also mask quiet sounds, making it harder to
understand them (Gelfand (2017); Lyon (2017)).

Therefore, the effects of reverberation on the HAS can vary widely, depending
on the duration and intensity of the phenomenon, the frequency content of the
sound source, the individual characteristics of the listener’s hearing, etc.
In general, it is desirable to optimize the amount of reverberation in a given
acoustic scenario to ensure the best possible listening experience for the
intended audience. In this situation, dereverberation is a process that can be
used to reduce or remove the effects of excessive reverberation from an audio
signal. This is typically used in situations where the audio signal has
already been recorded in a reverberant environment, and the resulting
reverberation is unwanted or detrimental to the quality of the recorded audio
(Naylor et al. (2010)).

In recent years, data-driven dereverberation methods have become increasingly
popular due to their ability to learn complex mappings between the anechoic
and reverberant signals. While showing promising results, these methods also
have a wide range of applicability, including speech recognition, speaker
verification, music production, etc (Xu et al. (2014); Hershey et al. (2016)).

Still, an important challenge for these methods is the need for large amounts
of training data, with paired clean and reverberant audio signals for
training, which, ideally, should cover a wide range of acoustic environments,
microphone types, and speaker characteristics, to ensure that a model will
generalize well to new scenarios. Examples of dereverberation datasets include
the REVERB challenge dataset (Kinoshita et al. (2013)), BUT Speech@FIT Reverb
Database (Szöke et al. (2019)), and VoiceBank-SLR (Fu et al. (2022)), among
others.

However, most of these datasets carry some limitations, such as narrow-band
anechoic data, i.e., only speech is considered as a signal of interest, and
absence of Impulse Response (IR) data, i.e., usually only anechoic and
reverberant data are paired for supervised training. Hence, in this paper, we
introduce a new data-compilation ensemble, primarily intended for training
data-driven dereverberation models capable of dealing with full bandwidth
audio signals, e.g., speech, song, music etc. We offer pairs of natural
anechoic and IR data, compiled from datasets licensed under Creative Commons
Attribution 4.0 International License, together with Python scripts for
convolution-based auralization, under the hypothesis that these ensembles
could serve as better training and evaluation tools for such algorithms.

The remainder of this paper is organized as follows: Section 2 details the
methodology used for selecting the anechoic and IR data, and then synthesizing
the auralized data. Section 3 describes the resulting auralized data and their
annotations. Finally, Section 4 presents an overall discussion of our obtained
results, together with some pertinent considerations to conclude our study.

##  2 Methods

To simulate or recreate an acoustic environment, such as a concert hall or a
recording studio, using computer algorithms and specialized software, one can
resort to a well known technique called auralization. It involves measurements
or simulations of the space’s physical properties to mimic its acoustic
characteristics. To proceed accordingly, there are various techniques
available, such as acoustic ray tracing and convolution-based auralization
(Kleiner et al. (1993)).

From a signal processing point of view, the convolution is a mathematical
operation that describes the interaction between two signals. In the context
of auralization, convolution can be used to simulate the effects of an
acoustic environment on an audio signal. The basic idea is to convolve an
audio signal with an IR that describes the acoustic characteristics of a room
or space. The resultant output produces a new signal that represents the
original audio signal after it has been modified by the acoustic
characteristics of the room (Allen and Berkley (1979); Oppenheim et al.
(1997)).

In the context of data-driven dereverberation methods, the choice for either
natural, synthetic or joint datasets, i.e., those obtained by processing
combinations of natural recordings with synthetic sounds, directly impacts in-
the-wild applicability. This occurs because, when models are trained on mostly
synthetic data, they usually don’t generalize well for real-world scenarios.
However, the majority of recent studies use joint combinations, e.g., by
convolving anechoic data with naturally recorded or synthesized IR data etc,
perhaps because it would be too cumbersome to naturally acquire all the
necessary data (dos Santos et al. (2022)). Hence our dataset attempts to reach
a balance between these features, compiling real anechoic and IR data, and
then synthetically producing auralized data ensembles by means of convolution.

###  2.1 Anechoic data

Trying to cover a wide range of full-bandwidth audio signals, we chose to
compile signals of interest from three different categories: speech, song, and
musical instruments.

Since RAVDESS (Livingstone and Russo (2018)) is a well-known dataset with
subsets of emotional speech and song, it was chosen to cover the first two
aforementioned categories. Its speech-only portion comprises $1,440$ samples
performed by $24$ actors ($12$ male and $12$ female), vocalizing $2$
lexically-matched statements ("kids are talking by the door" and "dogs are
sitting by the door") in a "neutral" North American accent. Each expression is
pronounced at $2$ different levels of emotional intensity (normal and strong),
with an additional neutral expression, resulting in a total amount of $60$
trials per actor. Speech emotions include calm, happy, sad, angry, fearful,
surprise, and disgust. Its song-only portion is quite similar and comprises
$1,012$ samples performed by $23$ actors, singing the same $2$ lexically-
matched statements. Song emotions include only neutral, calm, happy, sad,
angry, and fearful expressions. The original sample rate is fixed at 48 kHz,
and Figures 1 (a) and (b) illustrate histograms with the original duration of
files in each RAVDESS subset that we used.

Since the acoustic guitar is a popular musical instrument, for a variety of
reasons, including its ability to produce polyphonic sound and its musical
versatility, we chose to use a subset from GuitarSet (Xi et al. (2018)),
referred to as audio_mono-mic. It comprises $360$ samples performed by $6$
musicians, playing $30$ twelve to sixteen bar excerpts from lead-sheets in a
variety of keys, tempos, and musical genres. Recording was performed using a
Neumann U87 condenser microphone, placed at approximately 30 cm in front of
the 18th fret of the guitar. The original sample rate is fixed at 44.1 kHz,
and Figure 1 (c) illustrates a histogram with the original duration of files
in this GuitarSet subset.

####  2.1.1 Anechoic data processing

Considering the great difference between the duration of signals in RAVDESS
and GuitarSet, we chose to split the samples in GuitarSet into segments of
smaller duration, fixed at 5 s, resulting in $2,004$ different samples.
Another reason behind this decision is that when choosing the length of
anechoic data, it is important to strike a balance between the computational
cost of the convolution operation and the length of the segments. If they are
too short, the resulting audio signals may not capture the full extent of the
room’s acoustics, and if the segments are too long, the convolution operation
may become too computationally intensive.

Moreover, since the sample rates of RAVDESS and GuitarSet are different, we
also chose to up-sample the GuitarSet segments to 48 kHz, thus standardizing
this value for all files in our dataset.

###  2.2 IR data

IR data was curated from the Open Acoustic Impulse Response (Open AIR)
Library, which is an online database of Acoustic Impulse Response (AIR) data.
Since the original metadata in this database provides information about space
category, IR duration, etc., the samples were chosen in order to have a
balance between a selected variety of open and enclosed spaces, with IRs in
the range of a few milliseconds up to 5 s. Figure 1 (d) illustrates a
histogram with the original duration of the selected IR data. Altogether, $65$
IRs were chosen.

Figure 1: Original duration of files in the (a) speech, (b) song, (c) guitar,
and (d) IR subsets of our compilation.

####  2.2.1 IR data processing

The selected IR data was found in a variety of formats, e.g., B-format, MS
Stereo, mono, etc. And at sample rates varying from 44.1 kHz to 96 kHz. Since
all the anechoic data was compiled in mono format at 48 kHz, the selected IR
data was first converted to mono, then normalized to prevent files from
clipping, and finally, each IR was either down- or up-sampled to 48 kHz.

##  3 Results

The processed anechoic and IR data is hosted at Zenodo, with an approximate
total file size of 1.3 GB. For simplicity, all samples in our data-compilation
were renamed, e.g., $guitar\\_0000$, $rir\\_0000$, $song\\_0000$,
$speech\\_0000$, and so on. To synthesize different versions of the auralized
data ensemble, the reader is referred to GitHub, where Python scripts are
available for downloading the base data-compilation and synthesizing a chosen
version of the auralized data ensemble. Table 1 illustrates the differences
between all versions, detailing the number of song, speech, guitar, IR and
auralized samples in each one, together with their respective total file size
and duration.

Table 1: Number of anechoic, IR and resultant auralized data samples, together
with their respective total duration and file size for each ensemble version.

| Tiny | Small | Medium | Large
---|---|---|---|---
Song samples | 100 | 500 | 1,012 | 1,012
Speech samples | 100 | 500 | 1,012 | 1,440
Guitar samples | 100 | 500 | 1,012 | 2,004
IR samples | 5 | 9 | 33 | 65
Auralized samples | 1,500 | 13,500 | 100,188 | 289,640
Total duration | 3.2 h | 30.41 h | 221.77 h | 658.08 h
Total file size | 1.1 GB | 10.5 GB | 76.6 GB | 227.5 GB

###  3.1 Anechoic data annotation

Since the ancehoic data in our compilation is comprised of speech and song
samples with original annotations pertinent to emotion recognition, together
with acoustic guitar sounds with original annotations pertinent to Music
Information Retrieval (MIR), we provide metadata which can be used to trace
back each sample to its original annotations. This is done because we do not
intend for this dataset to be limited to dereverberation tasks only, but also
to be used for applications such as emotion recognition and MIR in more
challenging scenarios, i.e., in the presence of convolutive noise.

###  3.2 IR data annotation

Since the computational effort in dereverberation tasks is highly intertwined
with the RT values of IR data, $RT_{20}$ values were extracted from each IR
sample using ITA Toolbox. Figure 2 illustrates a histogram with the extracted
$RT_{20}$ values from the selected the IR data.

Figure 2: Extracted $RT_{20}$ values from the selected the IR data.

##  4 Discussion and Conclusion

Overall, the data-compilation ensemble presented in this work provides a
diverse and comprehensive set of acoustic scenes for use in dereverberation
tasks, as well as some other audio signal processing applications, such as
emotion recognition and MIR. By combining different types of signals of
interest, including speech, song, and acoustic guitar sounds, with a variety
of IRs, we provide a challenging dataset for researchers working on
dereverberation and related fields.

The dataset is available in different sizes, from a tiny version with limited
data, to a large version with almost $300,000$ samples, allowing users to
choose the most suitable version for their specific research needs. The
dataset also includes metadata that can be used to trace back each sample to
its original annotations, facilitating the use of the dataset for tasks such
as emotion recognition and MIR.

We hope that this dataset will be useful for researchers working on
dereverberation and related fields, and we encourage its use in future
research. We also believe that the diversity and variability of the dataset
can facilitate the development of more robust and generalizable algorithms for
dereverberation and other audio signal processing tasks.

## References

  * Beranek [2004] Leo Leroy Beranek.  _Concert halls and opera houses: music, acoustics, and architecture_ , volume 2.  Springer, 2004.
  * Lapointe [2010] Yannick Lapointe.  Understanding and crafting the mix: The art of recording.  _Intersections_ , 31(1):209, 2010.
  * Gelfand [2017] Stanley A Gelfand.  _Hearing: An introduction to psychological and physiological acoustics_.  CRC Press, 2017.
  * Lyon [2017] Richard F Lyon.  _Human and machine hearing: extracting meaning from sound_.  Cambridge University Press, 2017.
  * Naylor et al. [2010] Patrick A Naylor, Nikolay D Gaubitch, et al.  _Speech dereverberation_ , volume 2.  Springer, 2010.
  * Xu et al. [2014] Yong Xu, Jun Du, Li-Rong Dai, and Chin-Hui Lee.  A regression approach to speech enhancement based on deep neural networks.  _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , 23(1):7–19, 2014.
  * Hershey et al. [2016] John R Hershey, Zhuo Chen, Jonathan Le Roux, and Shinji Watanabe.  Deep clustering: Discriminative embeddings for segmentation and separation.  In _2016 IEEE international conference on acoustics, speech and signal processing (ICASSP)_ , pages 31–35. IEEE, 2016.
  * Kinoshita et al. [2013] Keisuke Kinoshita, Marc Delcroix, Takuya Yoshioka, Tomohiro Nakatani, Emanuel Habets, Reinhold Haeb-Umbach, Volker Leutnant, Armin Sehr, Walter Kellermann, Roland Maas, et al.  The reverb challenge: A common evaluation framework for dereverberation and recognition of reverberant speech.  In _2013 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics_ , pages 1–4. IEEE, 2013.
  * Szöke et al. [2019] Igor Szöke, Miroslav Skácel, Ladislav Mošner, Jakub Paliesek, and Jan Černockỳ.  Building and evaluation of a real room impulse response dataset.  _IEEE Journal of Selected Topics in Signal Processing_ , 13(4):863–876, 2019.
  * Fu et al. [2022] Szu-Wei Fu, Cheng Yu, Kuo-Hsuan Hung, Mirco Ravanelli, and Yu Tsao.  Metricgan-u: Unsupervised speech enhancement/dereverberation based only on noisy/reverberated speech.  In _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 7412–7416. IEEE, 2022.
  * Kleiner et al. [1993] Mendel Kleiner, Bengt-Inge Dalenbäck, and Peter Svensson.  Auralization-an overview.  _Journal of the Audio Engineering Society_ , 41(11):861–875, 1993.
  * Allen and Berkley [1979] Jont B Allen and David A Berkley.  Image method for efficiently simulating small-room acoustics.  _The Journal of the Acoustical Society of America_ , 65(4):943–950, 1979.
  * Oppenheim et al. [1997] Alan V Oppenheim, Alan S Willsky, Syed Hamid Nawab, and Jian-Jiun Ding.  _Signals and systems_ , volume 2.  Prentice hall Upper Saddle River, NJ, 1997.
  * dos Santos et al. [2022] Arthur dos Santos, Pedro de Oliveira, and Bruno Masiero.  A retrospective on multichannel speech and audio enhancement using machine and deep learning techniques.  In _Proceedings of the 24th International Congress on Acoustics_ , pages 173–184, Gyeongju, South Korea, October 2022.
  * Livingstone and Russo [2018] Steven R Livingstone and Frank A Russo.  The ryerson audio-visual database of emotional speech and song (ravdess): A dynamic, multimodal set of facial and vocal expressions in north american english.  _PloS one_ , 13(5):e0196391, 2018.
  * Xi et al. [2018] Qingyang Xi, Rachel M Bittner, Johan Pauwels, Xuzhou Ye, and Juan Pablo Bello.  Guitarset: A dataset for guitar transcription.  In _ISMIR_ , pages 453–460, 2018.