# Performance of Refined Isogeometric Analysis in Solving
Quadratic Eigenvalue Problems

Ali Hashemian<EMAIL_ADDRESS>Daniel Garcia  David Pardo  Victor M.
Calo  BCAM – Basque Center for Applied Mathematics, Bilbao, Basque Country,
Spain  IDAEA – Institute of Environmental Assessment and Water Research,
Barcelona, Cataluña, Spain  University of the Basque Country UPV/EHU, Leioa,
Basque Country, Spain  Ikerbasque – Basque Foundation for Sciences, Bilbao,
Basque Country, Spain  Faculty of Science and Engineering, Curtin University,
Perth, Australia

###### Abstract

Certain applications that analyze damping effects require the solution of
quadratic eigenvalue problems (QEPs). We use refined isogeometric analysis
(rIGA) to solve quadratic eigenproblems. rIGA discretization, while conserving
desirable properties of maximum-continuity isogeometric analysis (IGA),
reduces the interconnection between degrees of freedom by adding low-
continuity basis functions. This connectivity reduction in rIGA’s algebraic
system results in faster matrix LU factorizations when using multifrontal
direct solvers. We compare computational costs of rIGA versus those of IGA
when employing Krylov eigensolvers to solve quadratic eigenproblems arising in
2D vector-valued multifield problems. For large problem sizes, the
eigencomputation cost is governed by the cost of LU factorization, followed by
costs of several matrix–vector and vector–vector multiplications, which
correspond to Krylov projections. We minimize the computational cost by
introducing $C^{0}$ and $C^{1}$ separators at specific element interfaces for
our rIGA generalizations of the curl-conforming Nédélec and divergence-
conforming Raviart–Thomas finite elements. Let $p$ be the polynomial degree of
basis functions; the LU factorization is up to
${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$ times faster when using rIGA
compared to IGA in the asymptotic regime. Thus, rIGA theoretically improves
the total eigencomputation cost by ${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$
for sufficiently large problem sizes. Yet, in practical cases of moderate-size
eigenproblems, the improvement rate deteriorates as the number of computed
eigenvalues increases because of multiple matrix–vector and vector–vector
operations. Our numerical tests show that rIGA accelerates the solution of
quadratic eigensystems by ${\altmathcal{O}(p-1)}$ for moderately sized
problems when we seek to compute a reasonable number of eigenvalues.

###### keywords:

Quadratic eigenvalue problems; refined isogeometric analysis; computational
cost improvement; Krylov eigensolvers.

##  1 Introduction

Refined isogeometric analysis (rIGA), introduced by Garcia et al. [1], proved
to be a successful extension of isogeometric analysis (IGA) [2] when using
multifrontal direct solvers to approximate solutions to partial differential
equations (c.f., [3, 4, 5, 6, 7]). rIGA preserves the desirable properties of
maximum-continuity IGA discretizations while partitioning the computational
domain into macroelements interconnected by low-continuity basis functions.
This reduced connectivity exploits the recursive partitioning capability of
multifrontal direct solvers, significantly reducing the solution cost. In
practice, when using rIGA discretizations for gradient-conforming $H^{1}$
spaces, the matrix factorization is asymptotically ${\altmathcal{O}(p^{2})}$
times faster in large domains compared to IGA (where $p$ is the polynomial
order of B-spline bases). In this context, rIGA improves the performance of
discretizations by curl-conforming H(curl) and divergence-conforming H(div)
spaces by up to ${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$. In all cases, rIGA
also outperforms traditional finite element analysis (FEA) when considering a
fixed number of elements and using the same polynomial order (see [1, 3]).

IGA is a widely used numerical technique to solve eigenvalue problems (see,
e.g., [8, 9, 10, 11, 12, 13, 14, 15, 16, 17]). Recently, the application of
rIGA in eigenanalysis has been also investigated in [18, 6]. However, most of
the studied cases in the context of both maximum-continuity and refined IGA
frameworks are limited to the (linear) generalized eigenvalue problems, which
have the following form:

$\displaystyle(\textbf{A}-\lambda\textbf{B})\\!\>\textbf{u}=\textbf{0}\,,$
(1.1)

resulting in real eigenvalues when B is positive definite and A is Hermitian.
Quadratic eigenvalue problems (QEPs), commonly expressed as

$\displaystyle(\textbf{K}+\lambda\textbf{C}+\lambda^{2}\textbf{M})\\!\>\textbf{u}=\textbf{0}\,,$
(1.2)

are the most applicable subset of polynomial eigenvalue problems. They are of
great interest in multiple scientific and engineering case studies such as the
vibration of flexible mechanisms [19], fluid–structure interactions [20],
constrained least-squares [21], electromagnetic wave propagation in conductive
media [22], and acoustic fluid in a cavity with absorbing boundaries [23]. In
this work, we investigate the beneficial effect of using rIGA on
eigencomputation cost of solving QEPs. The arising quadratic eigensystems have
complex eigenvalues in many practical occasions. Thus, as opposed to the
generalized eigenproblems, some numerical algorithms (e.g., spectrum slicing,
see Section 4) are inapplicable to quadratic cases. Additionally, the number
of computed eigenvalues plays an important role in the overall
eigencomputation cost of QEPs. As a result, the implementation of rIGA
discretizations shows different levels of improvement in linear (generalized)
and quadratic eigenproblems (c.f., [6] for rIGA cost improvements in solving
generalized eigenproblems).

When solving a quadratic eigenproblem, we may compute the eigensolution by
linearizing the original eigensystem. Frequently used Krylov eigensolvers
commonly perform an Arnoldi recurrence to project the equivalent linearized
eigenproblem onto Krylov subspace111For some specific cases, quadratic
eigensolvers may alternatively employ the nonsymmetric or pseudo-Lanczos
algorithms (see, e.g, [24, 25]).. Then, they find a finite number of
eigenvalues using a Rayleigh–Ritz approximation (see, e.g., [26, 27, 24, 25]).
The projection process entails one (or more) LU factorization and multiple
matrix–vector operations — which are not naturally as expensive as matrix
factorization but should be taken into account because the eigencomputation
entails a significant number of these operations. Eigencomputations are among
the most expensive numerical approaches, especially for multidimensional
systems and moderate to large problem sizes when we seek to compute many
eigenvalues. Hashemian et al. [6] identify the three most expensive numerical
operations of solving generalized Hermitian eigenproblems, namely: matrix
factorization, forward/backward eliminations (i.e., multiplications of LU
factors by vectors), and matrix–vector products (i.e., multiplications of
system matrices by vectors). They show for eigenproblems in $H^{1}$ that, when
using multifrontal direct solvers, rIGA improves matrix factorization
asymptotically by ${\altmathcal{O}(p^{2})}$ for a sufficiently large number of
degrees of freedom. The improvement factor of forward/backward eliminations is
asymptotically up to $\altmathcal{O}(p)$. At the same time, the higher number
of nonzero entries in the rIGA matrices slightly degrades the matrix–vector
product cost. These improvement/degradation ratios also hold for quadratic
eigensystems. The difference, however, is in the number of times the
eigensolver calls each operation. In addition to the mentioned numerical
operations, the vector–vector product plays an important role in
eigencomputation cost of QEPs. Despite being considered as an inherently cheap
numerical operation, it has an important contribution to the total cost when
computing a large number of eigenvalues (see Section 5).

Herein, we use rIGA discretizations to solve quadratic eigensystems of
multidimensional problems with large numbers of degrees of freedom. We
consider two vector-valued multifield model problems arising in
electromagnetic wave propagation in coductive media, and vibration of an
acoustic fluid contained in a cavity with absorbing walls (see Section 2). We
use higher-order B-spline generalizations of the curl-conforming Nédélec and
divergence-conforming Raviart–Thomas finite elements (see, e.g., [10, 28, 29,
30]) to discretize our 2D electromagnetic and vibroacoustic problems,
respectively. These spaces allow us to obtain solution fields free of spurious
eigenmodes (i.e., modes with no physical meaning). We briefly describe IGA and
rIGA discretizations in Section 3. Section 4 describes the eigensolution
procedure applied to QEPs. In Section 5, we theoretically estimate the
eigencomputation costs based on the numerical operations employed by the
quadratic eigensolver. We show that rIGA improves the efficiency of the most
expensive numerical operations of the eigenanalysis. Our theoretical analysis
shows that an improvement of ${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$ in
eigencomputation cost is asymptotically possible when employing rIGA in
multifield problems discretized by H(curl) and H(div) spaces. Indeed, for
sufficiently large problems, the matrix factorization governs the solution
cost. However, in practical moderate-size problems, the numerical tests of
Section 6 show that rIGA reduces the eigencomputation cost by a factor of
approximately ${\altmathcal{O}(p-1)}$ when computing a reasonable number of
eigenvalues. Finally, we draw the main conclusions in Section 7.

##  2 Model problems

We consider the eigensolutions of two different quadratic eigenvalue problems.
The first eigenproblem arises in electromagnetic wave propagation; the second
one arises from the analysis of an acoustic (i.e., inviscid, compressible,
barotropic) fluid contained in a cavity with absorbing walls.

###  2.1 Electromagnetic wave propagation in a conductive medium

Let us define the H(curl)-conforming functional spaces on the computational
domain ${\Omega\subset\mathbb{R}^{2}}$ as

$\displaystyle\textit{{H}}({\rm curl};\Omega)$
$\displaystyle\coloneqq\bigg{\\{}\textbf{W}={\textstyle\footnotesize\begin{bmatrix}W_{x}\\\
W_{y}\end{bmatrix}}\in\big{(}L^{2}(\Omega)\big{)}^{2}:\nabla\times\textbf{W}\in\big{(}L^{2}(\Omega)\big{)}^{2}\bigg{\\}}\,,$
(2.1) $\displaystyle\textit{{H}}_{0}({\rm curl};\Omega)$
$\displaystyle\coloneqq\bigg{\\{}\textbf{W}\in\textit{{H}}({\rm
curl};\Omega):\textbf{W}\times\textbf{n}=\textbf{0}{\rm~{}on~{}}\partial\Omega\bigg{\\}}\,,$
(2.2)

where n is the outward unit normal vector on the boundary ${\partial\Omega}$.
We consider the eigenvalue problem of the electromagnetic wave propagation
equation in a conductive medium:

$\displaystyle\left\\{\begin{matrix}[l]{\rm
Find~{}}\omega\in\mathbb{C}{\rm~{}and~{}}\textbf{E}={\textstyle\footnotesize\begin{bmatrix}E_{x}\\\
E_{y}\end{bmatrix}}\in\textit{{H}}_{0}({\rm
curl};\Omega),{\rm~{}with~{}}\textbf{E}:\Omega\rightarrow\mathbb{C}^{2},{\rm~{}such~{}that}\\\\[7.0pt]
\qquad\begin{aligned}
\nabla\times(\,\tilde{\mu}^{-1}\nabla\times\textbf{E})+\tilde{\boldsymbol{\sigma}}\textbf{E}&=\textbf{0}\,,\quad&&{\rm
in~{}}\Omega\,,\\\\[3.0pt]
\textbf{E}\times\textbf{n}&=\textbf{0}\,,\quad&&{\rm
on~{}}\partial{\Omega}\,,\end{aligned}\end{matrix}\right.$ (2.3)

where $\omega$ is the eigenfrequency, E is the electromagnetic field,
${\tilde{\boldsymbol{\sigma}}\coloneqq\boldsymbol{\sigma}+{\rm
i\\!\>}\omega\varepsilon\\!\>\textbf{I}}$ and ${\tilde{\mu}\coloneqq{\rm
i\\!\>}\omega\mu}$ with $\boldsymbol{\sigma}$ being the conductivity matrix,
${\rm i\\!\>}$ being the imaginary unit, and $\varepsilon$ and $\mu$ being the
electric permittivity and magnetic permeability, respectively. We rewrite the
strong form of (2.3) as

$\displaystyle\nabla\times(-\mu^{-1}\nabla\times\textbf{E})-{\rm
i\\!\>}\omega\\!\>\boldsymbol{\sigma}\textbf{E}+\varepsilon\omega^{2}\textbf{E}=\textbf{0}\,,$
(2.4)

and consider ${\textbf{W}\in\textit{{H}}_{0}({\rm curl};\Omega)}$ as an
arbitrary test function to construct the sesquilinear forms

$\displaystyle a\\!\>(\textbf{W},\textbf{E})$
$\displaystyle\coloneqq\rmintop\nolimits_{\Omega}(\overline{\nabla\times\textbf{W}})\cdot(\nabla\times-\mu^{-1}\textbf{E})\,d\Omega\,,$
(2.5) $\displaystyle b\\!\>(\textbf{W},\textbf{E})$
$\displaystyle\coloneqq\rmintop\nolimits_{\Omega}\overline{\textbf{W}}\cdot\varepsilon\\!\>\textbf{E}\,d\Omega\,,$
(2.6) $\displaystyle c\\!\>(\textbf{W},\textbf{E})$
$\displaystyle\coloneqq\rmintop\nolimits_{\Omega}\overline{\textbf{W}}\cdot(-\\!\>{\rm
i\\!\>}\boldsymbol{\sigma}\textbf{E})\,d\Omega\,,$ (2.7)

where overbar means the conjugate form of the complex vector space and $\cdot$
denotes the inner product. We build the weak form of (2.3) as follows:

$\displaystyle
a\\!\>(\textbf{W},\textbf{E})+\omega\\!\>c\\!\>(\textbf{W},\textbf{E})+\omega^{2}b\\!\>(\textbf{W},\textbf{E})=0\,.$
(2.8)

Let us define

$\displaystyle\altmathcalb{V}_{0,h}^{\\!\>{\rm
curl}}(\Omega)\coloneqq\left\\{\textbf{W}_{h}\in\textit{{H}}({\rm
curl};\Omega)\cap\textit{{H}}^{1}(\Omega):\textbf{W}_{h}\times\textbf{n}=\textbf{0}{\rm~{}on~{}}\partial{\Omega}\right\\},$
(2.9)

as our curl-conforming approximate space, where index $h$ refers to
approximate values obtained by discretization. We use a Galerkin
discretization of the continuous eigenproblem and express the approximate
solution field as follows:

$\displaystyle\textbf{E}_{h}=\rmsum_{i=0}^{N-1}\boldsymbol{\phi}_{i}^{\\!\>{\rm
curl}}\tilde{E}_{h,i}\,,$ (2.10)

where $\tilde{E}_{h,i}$ represents the ${i\\!\>}$-th nodal component of the
discrete solution $\tilde{\textbf{E}}_{h}$ given by

$\displaystyle\tilde{\textbf{E}}_{h}=\big{[}\tilde{E}_{x,0}\\!\>,...\\!\>,\tilde{E}_{x,N_{x}-1},\tilde{E}_{y,0}\\!\>,...\\!\>,\tilde{E}_{y,N_{y}-1}\big{]}\,,$
(2.11)

and ${\boldsymbol{\phi}_{i}^{\\!\>{\rm curl}}}$ is its corresponding curl-
conforming vector-valued basis function. Herein, ${N\coloneqq N_{x}+N_{y}}$ is
the total number of degrees of freedom, being $N_{x}$ and $N_{y}$ the numbers
of degrees of freedom of vector fields in $x$ and $y$ directions,
respectively. Thus, we write the approximate form of (2.8) as

$\displaystyle
a\\!\>(\textbf{W}_{h}\\!\>,\textbf{E}_{h})+\omega_{h}\\!\>c\\!\>(\textbf{W}_{h}\\!\>,\textbf{E}_{h})+\omega_{h}^{2}\\!\>b\\!\>(\textbf{W}_{h}\\!\>,\textbf{E}_{h})=0\,,$
(2.12)

that results in the following quadratic eigenvalue problem in the discrete
form:

$\displaystyle\big{(}\textbf{K}+\omega_{h}\textbf{C}+\omega_{h}^{2}\\!\>\textbf{M}\big{)}\\!\>\tilde{\textbf{E}}_{h}=\textbf{0}\,.$
(2.13)

###### Remark 1

In the context of eigenvalue analysis, we refer to the discrete solution
vector $\tilde{\textbf{E}}_{h}$ as the eigenvector and the continuous solution
filed (either exact E or approximate $\textbf{E}_{h}$) as the eigenfunction.
For an $N$-degree-of-freedom system, the eigenproblem (2.13) has $2N$
approximate eigenpairs referring to either ${(\omega_{h},\textbf{E}_{h})}$ or
${(\omega_{h},\tilde{\textbf{E}}_{h})}$ .

The stiffness and mass matrices in eigenproblem (2.13), K and M, respectively,
are symmetric and M is positive definite. The damping matrix C is skew-
Hermitian (i.e., ${\textbf{C}=-\textbf{C}^{T}}$, where superscript $T$ denotes
the conjugate transpose). Thus, the current model problem results in a
gyroscopic eigensystem (c.f., [31]) with a Hamiltonian spectrum structure
whose eigenvalues are symmetric with respect to the imaginary axis. We refer
to ${{\rm Re}\\!\>(\omega_{h})}$ and ${{\rm Im}\\!\>(\omega_{h})\geq 0}$ as
the angular frequency and exponential decay of the electromagnetic wave with
time, respectively [22].

###### Remark 2

In a non-conductive medium with ${\boldsymbol{\sigma}=\textbf{0}}$ , the
eigensystem (2.13) is reduced to the generalized eigenproblem
${(\textbf{A}-\lambda_{h}\textbf{B})\\!\>\tilde{\textbf{E}}_{h}=\textbf{0}}$
and commonly referred to as the Maxwell’s eigenvalue problem (i.e., the curl-
curl operator eigenproblem [10, 32]). We set
${\lambda_{h}\coloneqq\omega_{h}^{2}}$ and define
${\textbf{A}\coloneqq-\textbf{K}}$ and ${\textbf{B}\coloneqq\textbf{M}}$ to
follow the notations of generalized eigenproblems described in Section 4. In
this case, since B is positive definite and A is positive semidefinite, all
$\lambda_{h}$ values are real non-negative and the electromagnetic wave does
not decay.

###  2.2 Acoustic fluid contained in a cavity with absorbing walls

We consider an acoustic fluid contained in a rigid cavity with some (or all)
of its walls covered by a viscoelastic material absorbing the acoustic energy
of the fluid [23]. Let us consider the domain ${\Omega\subset\mathbb{R}^{2}}$
and its boundary ${\partial\Omega\coloneqq\Gamma_{A}\cup\Gamma_{R}}$ , being
$\Gamma_{A}$ and $\Gamma_{R}$ the union of absorbing and rigid boundaries,
respectively. We define the H(div)-conforming functional spaces as follows:

$\displaystyle\textit{{H}}({\rm div};\Omega)$
$\displaystyle\coloneqq\bigg{\\{}\textbf{w}={\textstyle\footnotesize\begin{bmatrix}w_{x}\\\
w_{y}\end{bmatrix}}\in\big{(}L^{2}(\Omega)\big{)}^{2}:\nabla\cdot\textbf{w}\in
L^{2}(\Omega)\bigg{\\}}\,,$ (2.14) $\displaystyle\textit{{H}}_{0}({\rm
div};\Omega)$ $\displaystyle\coloneqq\bigg{\\{}\textbf{w}\in\textit{{H}}({\rm
div};\Omega):\textbf{w}\cdot\textbf{n}\in
L^{2}(\partial\Omega){\rm~{}and~{}}\textbf{w}\cdot\textbf{n}=0~{}{\rm
on}~{}\Gamma_{R}\bigg{\\}}\,.$ (2.15)

We find the damped vibration modes of the fluid as complex solutions to the
following problem:

$\displaystyle\left\\{\begin{matrix}[l]{\rm
Find~{}}\lambda\in\mathbb{C},\textbf{u}={\textstyle\footnotesize\begin{bmatrix}u_{x}\\\
u_{y}\end{bmatrix}}\in\textit{{H}}_{0}({\rm div};\Omega){\rm~{}and~{}}P\in
H_{0}^{1}(\Omega),{\rm~{}with~{}}\textbf{u}:\Omega\rightarrow\mathbb{C}^{2}{\rm~{}and~{}}P:\Omega\rightarrow\mathbb{C},{\rm~{}such~{}that}\\\\[7.0pt]
\qquad\begin{aligned} \rho\lambda^{2}\\!\>\textbf{u}+\nabla P&=0\,,\quad&&{\rm
in~{}}\Omega\,,\\\\[2.0pt] -\rho
c^{2}\\!\>\nabla\cdot\textbf{u}&=P\,,\quad&&{\rm in~{}}\Omega\,,\\\\[2.0pt]
(\alpha+\lambda\\!\>\beta)\\!\>\textbf{u}\cdot\textbf{n}&=P\,,\quad&&{\rm
on~{}}\Gamma_{A}\,,\\\\[2.0pt] \textbf{u}\cdot\textbf{n}&=0\,,\quad&&{\rm
on~{}}\Gamma_{R}\,,\end{aligned}\end{matrix}\right.$ (2.16)

where $P$ is the fluid pressure, u is the displacement field, $\rho$ is the
fluid density, and $c$ is the acoustic speed. Herein, ${{\rm
Im}\\!\>(\lambda)}$ and ${{\rm Re}\\!\>(\lambda)}$ are the vibration’s angular
frequency and decay rate, respectively, while coefficients $\alpha$ and
$\beta$ are related to the impedance of the viscoelastic material.

By considering an arbitrary test function ${\textbf{w}\in\textit{{H}}_{0}({\rm
div};\Omega)}$ and defining the sesquilinear forms

$\displaystyle\breve{a}\\!\>(\textbf{w},\textbf{u})$
$\displaystyle\coloneqq\rmintop\nolimits_{\Omega}(\nabla\cdot\overline{\textbf{w}})\\!\>(\,\rho
c^{2}\nabla\cdot\textbf{u})\,d\Omega+\rmintop\nolimits_{\Gamma_{A}}(\overline{\textbf{w}}\cdot\textbf{n})\\!\>(\alpha\\!\>\textbf{u}\cdot\textbf{n})\,d\\!\>\Gamma_{A}\,,$
(2.17) $\displaystyle\breve{b}\\!\>(\textbf{w},\textbf{u})$
$\displaystyle\coloneqq\rmintop\nolimits_{\Omega}\overline{\textbf{w}}\cdot\rho\\!\>\textbf{u}\,d\Omega\,,$
(2.18) $\displaystyle\breve{c}\\!\>(\textbf{w},\textbf{u})$
$\displaystyle\coloneqq\rmintop\nolimits_{\Gamma_{A}}(\overline{\textbf{w}}\cdot\textbf{n})\\!\>(\,\beta\\!\>\textbf{u}\cdot\textbf{n})\,d\\!\>\Gamma_{A}\,,$
(2.19)

we obtain the weak form of (2.16) involving only displacement variables as
follows (we use the ${\,\breve{}\,}$ symbol to distinguish between
sesquilinear forms of our two model problems):

$\displaystyle\breve{a}\\!\>(\textbf{w},\textbf{u})+\lambda\\!\>\breve{c}\\!\>(\textbf{w},\textbf{u})+\lambda^{2}\\!\>\breve{b}\\!\>(\textbf{w},\textbf{u})=0\,.$
(2.20)

We write another variational formulation in terms of the pressure that leads
to a nonlinear eigenproblem [23] that is not studied here. Let us consider

$\displaystyle\altmathcalb{V}_{0,h}^{\\!\>{\rm
div}}(\Omega)\coloneqq\left\\{\textbf{w}_{h}\in\textit{{H}}({\rm
div};\Omega)\cap\textit{{H}}^{1}(\Omega):\textbf{w}_{h}\cdot\textbf{n}\in
L^{2}(\partial\Omega){\rm~{}and~{}}\textbf{w}_{h}\cdot\textbf{n}=0{\rm~{}on~{}}\Gamma_{R}\right\\},$
(2.21)

as our divergence-conforming approximate space. We discretize the continuous
displacement field as

$\displaystyle\textbf{u}_{h}=\rmsum_{i=0}^{N-1}\boldsymbol{\phi}_{i}^{\\!\>{\rm
div}}\\!\>\tilde{u}_{h,i}\,,$ (2.22)

where ${\boldsymbol{\phi}_{i}^{\\!\>{\rm div}}}$ is the divergence-conforming
vector-valued basis and $\tilde{u}_{h,i}$ is its corresponding component of
the discrete solution field
${\tilde{\textbf{u}}_{h}=\big{[}\tilde{u}_{x,0}\\!\>,...\\!\>,\tilde{u}_{x,N_{x}-1},\tilde{u}_{y,0}\\!\>,...\\!\>,\tilde{u}_{y,N_{y}-1}\big{]}}$.
Thus, we write the approximate form of (2.20) as

$\displaystyle\breve{a}\\!\>(\textbf{w}_{h}\\!\>,\textbf{u}_{h})+\lambda_{h}\\!\>\breve{c}\\!\>(\textbf{w}_{h}\\!\>,\textbf{u}_{h})+\lambda_{h}^{2}\\!\>\breve{b}\\!\>(\textbf{w}_{h}\\!\>,\textbf{u}_{h})=0\,,$
(2.23)

resulting in

$\displaystyle\big{(}\breve{\textbf{K}}+\lambda_{h}\breve{\textbf{C}}+\lambda_{h}^{2}\\!\>\breve{\textbf{M}}\big{)}\\!\>\tilde{\textbf{u}}_{h}=\textbf{0}\,,$
(2.24)

which is a Hermitian quadratic eigenproblem (all system matrices are real
symmetric). The mass matrix $\breve{\textbf{M}}$ is positive definite, while
the stiffness $\breve{\textbf{K}}$ and damping $\breve{\textbf{C}}$ are
positive semidefinite matrices resulting in a stable vibration response, i.e.,
${{\rm Re}\\!\>(\lambda_{h})\leq 0}$. We compute the eigenvalues as a
combination of purely real and conjugate complex pairs. The former corresponds
to an overdamped non-oscillatory response, while the latter amounts to an
underdamped vibration.

##  3 Refined isogeometric analysis

We now consider multifield vectorial solutions and discretize the continuous
eigenproblems (2.3) and (2.16) using B-spline generalizations of curl-
conforming Nédélec and divergence-conforming Raviart–Thomas spaces,
respectively (see, e.g., [3]). We herein review some basic concepts of
maximum-continuity and refined IGA discretizations.

###  3.1 Bivariate B-spline space

Given the parametric domain
${\hat{\Omega}:\xi,\eta\in(0,1)^{2}\subset\mathbb{R}^{2}}$, we introduce the
bivariate B-spline space ${\altmathcal{S}^{p_{x},p_{y}}_{k_{x},k_{y}}}$ as

$\displaystyle\altmathcal{S}^{p_{x},p_{y}}_{k_{x},k_{y}}\coloneqq{\rm
span}\left\\{B_{i,j}^{\\!\>p_{x},p_{y}}\right\\}_{i=0,j=0}^{n_{x}-1,n_{y}-1},$
(3.1)

where $n$, $p$, and $k$ with their indices are the number of degrees of
freedom, polynomial degree, and continuity of basis functions in each
direction, respectively. The bivariate basis functions
${B_{i,j}^{\\!\>p_{x},p_{y}}}$ are given in the tensor product sense by:

$\displaystyle B_{i,j}^{\\!\>p_{x},p_{y}}\coloneqq
B_{i}^{\\!\>p_{x}}(\xi)\\!\>B_{j}^{\\!\>p_{y}}(\eta)\,,\quad
i=0,1,...,n_{x}-1\,,\quad j=0,1,...,n_{y}-1\,,$ (3.2)

where the univariate bases span their respective knot sequences

$\displaystyle\Xi$
$\displaystyle=[\underbrace{0,0,...,0}_{p_{x}+1},\xi_{p_{x}+1},\xi_{p_{x}+2}...,\xi_{n_{x}-1},\underbrace{1,1,...,1}_{p_{x}+1}]\,,$
(3.3) $\displaystyle\rm{H}$
$\displaystyle=[\underbrace{0,0,...,0}_{p_{y}+1},\eta_{p_{y}+1},\eta_{p_{y}+2}...,\eta_{n_{y}-1},\underbrace{1,1,...,1}_{p_{y}+1}]\,.$
(3.4)

We use the Cox–De Boor recursion formula [33] to obtain the univariate basis
functions, which in terms of $\xi$ reads:

$\displaystyle B_{i}^{\\!\>0}(\xi)$
$\displaystyle=\begin{cases}\begin{aligned}
&1\,,&&\xi_{i}\leq\xi<\xi_{i+1}\,,\\\ &0\,,&&{\rm
otherwise}\,,\end{aligned}\end{cases}$ (3.5) $\displaystyle
B_{i}^{\\!\>p_{x}}(\xi)$
$\displaystyle=\dfrac{\xi-\xi_{i}}{\xi_{i+p_{x}}-\xi_{i}}B_{i}^{\\!\>p_{x}-1}(\xi)+\dfrac{\xi_{i+p_{x}+1}-\xi}{\xi_{i+p_{x}+1}-\xi_{i+1}}B_{i+1}^{\\!\>p_{x}-1}(\xi)\,.$
(3.6)

When evaluating basis functions, we find the corresponding nonzero knot span
in (LABEL:eq.BasisFunBi0) and efficiently evaluate (3.6) avoiding any division
by zero. More details are given in [33, Algorithms A2.1 and A2.2].

###  3.2 IGA discretization

In the IGA context, we consider single multiplicities for all interior knots
in (3.3) and (3.4), providing maximum continuity of basis functions. Let us
consider polynomial degree $p$ for discretization and assume the 2D domain has
an ${n_{e}\times n_{e}}$ grid, being $n_{e}$ the number of elements in each
direction. We construct the curl- and divergence-conforming discrete spaces in
the parametric domain (c.f., [10, 3]) as

$\displaystyle\hat{\altmathcalb{V}}_{h,{\rm IGA}}^{\,{\rm
curl}}(\hat{\Omega})$
$\displaystyle\coloneqq\altmathcal{S}^{p-1,p}_{k-1,k}\times\altmathcal{S}^{p,p-1}_{k,k-1}\,,$
(3.7) $\displaystyle\hat{\altmathcalb{V}}_{h,{\rm IGA}}^{\,{\rm
div}}(\hat{\Omega})$
$\displaystyle\coloneqq\altmathcal{S}^{p,p-1}_{k,k-1}\times\altmathcal{S}^{p-1,p}_{k-1,k}\,,$
(3.8)

where ${k=p-1}$. The following vector-valued bivariate B-spline bases
characterize these spaces:

$\displaystyle\hat{\boldsymbol{\phi}}^{\\!\>{\rm curl}}_{\\!\>{\rm IGA}}$
$\displaystyle\coloneqq\begin{bmatrix}B_{0,0}^{\\!\>p-1,p}&...&B_{n_{e}+p-2,\\!\>n_{e}+p-1}^{\\!\>p-1,p}&0&...&0\\\\[3.0pt]
0&...&0&B_{0,0}^{\\!\>p,p-1}&...&B_{n_{e}+p-1,\\!\>n_{e}+p-2}^{\\!\>p,p-1}\end{bmatrix}\\!\>,$
(3.9) $\displaystyle\hat{\boldsymbol{\phi}}^{\\!\>{\rm div}}_{\\!\>{\rm IGA}}$
$\displaystyle\coloneqq\begin{bmatrix}B_{0,0}^{\\!\>p,p-1}&...&B_{n_{e}+p-1,\\!\>n_{e}+p-2}^{\\!\>p,p-1}&0&...&0\\\\[3.0pt]
0&...&0&B_{0,0}^{\\!\>p-1,p}&...&B_{n_{e}+p-2,\\!\>n_{e}+p-1}^{\\!\>p-1,p}\end{bmatrix}\\!\>.$
(3.10)

Fig. 1 illustrates examples of H(curl) and H(div) IGA discrete spaces as well
as univariate basis functions of the respective vector fields in the
parametric domain. For 2D problems, H(curl) and H(div) are the same spaces
rotated by 90 degrees in the parametric domain. Herein, for the sake of
simplicity, we use u in reference to the solution fields of both model
problems described in Section 2.

(a) H(curl) IGA discrete space

(b) H(div) IGA discrete space

Figure 1: Examples of H(curl) and H(div) spaces discretized by maximum-
continuity IGA with uniform ${16\times 16}$ elements, polynomial degree
${p=4}$, and continuity ${k=3}$. We show the univariate bases of the solution
fields $\hat{u}_{h,x}$ and $\hat{u}_{h,y}$ in the parametric domain, where the
blue and green color codes refer to basis functions of degree $p$ and ${p-1}$,
respectively.

Let us define ${\textbf{F}:\hat{\Omega}\rightarrow\Omega}$ as the geometric
mapping from the parametric space onto the physical domain and $D\textbf{F}$
as its Jacobian. We introduce our discrete spaces in the physical domain as

$\displaystyle\altmathcalb{V}_{h,{\rm IGA}}^{\,{\rm curl}}(\Omega)$
$\displaystyle\coloneqq\left\\{\textbf{u}_{h}\in\textit{{H}}({\rm
curl};\Omega)\cap\textit{{H}}^{1}(\Omega):\iota^{\\!\>{\rm
curl}}(\textbf{u}_{h})=\hat{\textbf{u}}_{h}\in\hat{\altmathcalb{V}}_{h,{\rm
IGA}}^{\,{\rm curl}}(\hat{\Omega})\right\\}\,,$ (3.11)
$\displaystyle\altmathcalb{V}_{h,{\rm IGA}}^{\,{\rm div}}(\Omega)$
$\displaystyle\coloneqq\left\\{\textbf{u}_{h}\in\textit{{H}}({\rm
div};\Omega)\cap\textit{{H}}^{1}(\Omega):\iota^{\\!\>{\rm
div}}(\textbf{u}_{h})=\hat{\textbf{u}}_{h}\in\hat{\altmathcalb{V}}_{h,{\rm
IGA}}^{\,{\rm div}}(\hat{\Omega})\right\\}\,,$ (3.12)

considering the following curl- and divergence-preserving mappings [3, 10]:

$\displaystyle\iota^{\\!\>{\rm curl}}(\textbf{u}_{h})$
$\displaystyle\coloneqq(D\textbf{F})^{T}(\textbf{u}_{h}\circ\textbf{F})\,,$
(3.13) $\displaystyle\iota^{\\!\>{\rm div}}(\textbf{u}_{h})$
$\displaystyle\coloneqq\det\\!\>(D\textbf{F})\\!\>(D\textbf{F})^{-1}(\textbf{u}_{h}\circ\textbf{F})\,.$
(3.14)

###  3.3 rIGA discretization

Refined isogeometric analysis (rIGA) is a discretization technique that
optimizes the performance of direct solvers. In particular, rIGA preserves the
optimal convergence order of direct solvers for a fixed number of elements in
the domain. Garcia et al. [1] first presented this strategy for $H^{1}$ spaces
and then extended it to H(curl), H(div), and $L^{2}$ spaces [3]. Starting from
the maximum-continuity IGA discretization, rIGA reduces the continuity of
certain basis functions by increasing the multiplicity of the respective
existing knots. Hence, we subdivide the computational domain into high-
continuity macroelements interconnected by low-continuity hyperplanes. These
hyperplanes coincide with the locations of the separators at different
partitioning levels of the multifrontal direct solvers. Thus, when solving
systems of equations arising in $H^{1}$-discretized problems, rIGA performs
matrix factorization asymptotically ${\altmathcal{O}(p^{2})}$ times faster in
large domains — and $\altmathcal{O}(p)$ faster in small domains — compared
to IGA. The improvement factors are ${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$
and ${\altmathcal{O}(p-1)}$ for large and small problems, respectively, when
considering H(curl) and H(div) spaces (mainly because we reduce to degree
${p-1}$ at certain spatial directions). In comparison to traditional FEA with
the same number of elements, rIGA provides even larger improvements.
Additionally, rIGA reduces the forward/backward elimination cost
asymptotically by up to $\altmathcal{O}(p)$ in $H^{1}$ and
${\altmathcal{O}(p-1)}$ in H(curl) and H(div) for sufficiently large domains
since the LU (or Cholesky) factors have fewer nonzero terms for rIGA. Consider
a mesh with a fixed element number in each direction being a power of two;
these improvements correspond to macroelements with 16 elements in each
direction for moderately sized problems (see [1, 3, 6]).

We discretize our multifield problems using H(curl) and H(div) spaces; thus,
we preserve the commutativity of the de Rham diagram [34] by reducing the
continuity in ${k-1}$ degrees. We achieve the commutativity using both
$C^{\\!\>0}$ and $C^{1}$ hyperplanes and reduce the continuity across
interfaces between subdomains (i.e., macroelements). Thus, we represent the
rIGA spaces in the parametric domain as follows:

$\displaystyle\hat{\altmathcalb{V}}_{h,{\rm rIGA}}^{\,{\rm
curl}}(\hat{\Omega})$
$\displaystyle\coloneqq\altmathcal{S}^{p-1,p}_{(k-1,0\\!\>|_{\\!\>\rm
vs}),\\!\>(k,1|_{\\!\>\rm hs})}\times\altmathcal{S}^{p,p-1}_{(k,1|_{\\!\>\rm
vs}),\\!\>(k-1,0\\!\>|_{\\!\>\rm hs})}\,,$ (3.15)
$\displaystyle\hat{\altmathcalb{V}}_{h,{\rm rIGA}}^{\,{\rm
div}}(\hat{\Omega})$
$\displaystyle\coloneqq\altmathcal{S}^{p,p-1}_{(k,1|_{\\!\>\rm
vs}),\\!\>(k-1,0\\!\>|_{\\!\>\rm
hs})}\times\altmathcal{S}^{p-1,p}_{(k-1,0\\!\>|_{\\!\>\rm
vs}),\\!\>(k,1|_{\\!\>\rm hs})}\,,$ (3.16)

where “vs” and “hs” denote the vertical and horizontal separators,
respectively. Accordingly, one updates (3.9)–(3.14) in a similar manner to
write rIGA basis functions, mappings, and discrete spaces in the physical
domain. Fig. 2 depicts an example of an rIGA discretization of the discrete
spaces of Fig. 1 obtained by two levels of symmetric partitioning, which
results in ${4\times 4}$ macroelements.

(a) H(curl) rIGA discrete space

(b) H(div) rIGA discrete space

Figure 2: H(curl) and H(div) rIGA spaces of the ${16\times 16}$ domain of Fig.
1 with ${p=4}$ and ${k=3}$, after two levels of partitioning that generate
four blocks in each direction. rIGA reduces the continuity of basis functions
by ${k-1}$ degrees across the macroelement separators (we show the low-
continuity bases corresponding to the first and second partitioning levels in
dark red and black, respectively). Thin lines in the mesh skeleton denote the
high-continuity element interfaces, while thick lines illustrate the
macroelement boundaries.

\begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt 206.28247pt
0.0pt,clip]{Figs/HcurlN64p3L0Matrix}\end{overpic} (a) IGA matrix: 9,112
$\times$ 9,112

\begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt 206.28247pt
0.0pt,clip]{Figs/HcurlN64p3L0Reorder}\end{overpic} (b) IGA reordered matrix:
$N_{nz}=\,$1,090,240

\begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt 206.28247pt
0.0pt,clip]{Figs/HcurlN64p3L0Lowertril}\end{overpic} (c) IGA lower factor:
$N_{nz}=\,$3,142,952

\begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt 206.28247pt
0.0pt,clip]{Figs/HcurlN64p3L2Matrix}\end{overpic} (d) rIGA matrix: 10,804
$\times$ 10,804

\begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt 206.28247pt
0.0pt,clip]{Figs/HcurlN64p3L2Reorder}\end{overpic} (e) rIGA reordered matrix:
$N_{nz}=\,$1,217,968

\begin{overpic}[width=433.62pt,trim=0.0pt 0.0pt 206.28247pt
0.0pt,clip]{Figs/HcurlN64p3L2Lowertril}\end{overpic} (f) rIGA lower factor:
$N_{nz}=\,$2,117,582

Figure 3: Matrix patterns for a ${64\times 64}$ domain discretization using
IGA and rIGA H(curl) spaces with ${p=4\\!\>}$: (a,d) original system matrices,
(b,e) reordered matrices obtained by nested dissection permutation, and (c,f)
lower triangular Cholesky factors. We obtain an rIGA discretization by
performing two partitioning levels, creating ${16\times 16}$ macroelements.

We investigate how rIGA modifies the system matrices and the respective LU (or
Cholesky) factors; thus, we consider an H(curl) space discretized by
${64\times 64}$ elements and polynomial order ${p=4}$. Figs. 3(a) and 3(d)
show the IGA and rIGA matrix patterns, respectively. We use two levels of
partitioning that creates ${16\times 16}$ macroelements. We use the nested
dissection permutation [35] to obtain the reordered matrices (Figs. 3(b) and
3(e)). An rIGA-discretized system has more number of degrees of freedom and,
hence, a higher number of nonzeros ($N_{nz}$). However, the lower continuity
of the mesh at hyperplanes leads to lower interconnections in the matrix
connectivity graph. Thus, the factorization is obtained faster and the
respective factors have a lower number of nonzero terms (Figs. 3(c) and 3(f)
).

##  4 Numerical solution to quadratic eigenvalue problems

We review the numerical solution approach to quadratic eigenvalue problems.
For the sake of brevity and to follow the same notation for both model
problems introduced in Section 2, we denote the stiffness, damping, and mass
matrices by K, C, and M, respectively. We also drop the subscript $h$ and
consider ${(\lambda,\tilde{\textbf{u}})}$ as the approximate eigenvalue and
eigenvector pair. We review the generalized eigenproblems since we first
linearize the quadratic eigenvalue problem and then solve a generalized
eigensystem (see, e.g., [20, 23, 27, 24]).

###  4.1 Generalized eigenproblems

Let us consider the following generalized eigenproblem,

$\displaystyle(\textbf{A}-\lambda\textbf{B})\\!\>\tilde{\textbf{u}}=\textbf{0}\,,$
(4.1)

characterized by a linear matrix pencil
${\altmathcalb{L}(\lambda)\coloneqq\textbf{A}-\lambda\textbf{B}\\!\>}$. The
term generalized distinguishes (4.1) from the standard eigenproblem for which
B is the unit matrix. If both A and B are Hermitian matrices and B is
definite, $\altmathcalb{L}(\lambda)$ is a Hermitian definite pencil, and all
eigenvalues are real; otherwise, some eigenvalues may become complex. When
solving QEPs, via linearization, the linear pencil, on most occasions, is non-
Hermitian or indefinite, resulting in complex eigenvalues. For linear
Hermitian definite pencils, most eigensolvers obtain a portion of the
eigenspectrum by executing a sequence of numerical algorithms, among which
are: shift-and-invert spectral transformation, Arnoldi (or Lanczos)
decomposition, restarting, spectrum slicing, and deflation (see, e.g., [6]).
However, not all algorithms apply to the quadratic case as we may deal with
non-Hermitian (or indefinite) linearized pencils. In the following, we present
a short review of those algorithms essential to the efficiency study of the
eigencomputation of QEPs.

####  4.1.1 Shift-and-invert spectral transformation

When solving eigenproblem (4.1), an efficient way of computing eigenpairs is
to shift the spectrum by ${s\in\mathbb{C}}$ and solve the shift-and-invert
spectral transformed eigenproblem222We first transform (4.1) to
${(\textbf{A}-s\\!\>\textbf{B})\\!\>\tilde{\textbf{u}}=(\lambda-s)\\!\>\textbf{B}\tilde{\textbf{u}}}$
and, then, invert it to reach (4.2). (see, e.g., [36, 37, 38, 39, 26]),

$\displaystyle(\textbf{A}-s\\!\>\textbf{B})^{-1}\textbf{B}\tilde{\textbf{u}}=\theta\\!\>\tilde{\textbf{u}}\,,$
(4.2)

where ${\theta=1/(\lambda-s)}$. In large multidimensional eigensystems with
multiple clustered eigenvalues, this approach is fast and accurate when
calculating eigenvalues near the shift $s$. By defining the operator matrix,

$\displaystyle\textbf{S}\coloneqq(\textbf{A}-s\\!\>\textbf{B})^{-1}\textbf{B}=\altmathcalb{L}(s)^{-1}\textbf{B}\,,$
(4.3)

we express the shift-and-invert eigenproblem as
${\textbf{S}\\!\>\tilde{\textbf{u}}=\theta\\!\>\tilde{\textbf{u}}}$.

####  4.1.2 Arnoldi decomposition

To solve the eigenproblem (4.2), when B is positive definite, we employ the
Arnoldi decomposition, which can deal with either Hermitian or non-Hermitian
A. For a given shift $s$ and either Cholesky or LU factorization of
$\altmathcalb{L}(s)$, the $m$-step Arnoldi decomposition,

$\displaystyle\textbf{S}\\!\>\textbf{V}_{m}=\textbf{V}_{m}\textbf{H}_{m}+\beta\\!\>\textbf{v}_{m+1}\textbf{e}_{m}^{T}\,,$
(4.4)

projects the eigensystem (4.2) of size $N$ onto Krylov subspace of size $m$
${(m\ll N)}$ by reducing the ${N\times N}$ operator matrix S to an ${m\times
m}$ upper Hessenberg matrix $\textbf{H}_{m}$. Herein, $\textbf{e}_{m}$ is the
${m\\!\>}$th coordinate vector and the term
${\beta\\!\>\textbf{v}_{m+1}\textbf{e}_{m}^{T}}$ is the residual of the
$m$-step Arnoldi decomposition (see, e.g., [26, 40]). In the above equation,

$\displaystyle\textbf{H}_{m}\coloneqq\begin{bmatrix}h_{1,1}&h_{1,2}&\cdots&h_{1,m-2}&h_{1,m-1}&h_{1,m}\\\
h_{2,1}&h_{2,2}&\cdots&h_{2,m-2}&h_{2,m-1}&h_{2,m}\\\
&h_{3,2}&\cdots&h_{3,m-2}&h_{3,m-1}&h_{3,m}\\\
&&\ddots&\vdots&\vdots&\vdots\\\ &&&h_{m-1,m-2}&h_{m-1,m-1}&h_{m-1,m}\\\
&&&&h_{m,m-1}&h_{m,m}\end{bmatrix},$ (4.5)

and
${\textbf{V}_{m}\coloneqq[\textbf{v}_{1},\textbf{v}_{2},...,\textbf{v}_{m}]}$
is the ${N\times m}$ matrix of Arnoldi basis vectors. We obtain the components
of $\textbf{H}_{m}$ and Arnoldi bases using the recurrence given by Algorithm
1.

Algorithm 1 Arnoldi decomposition for solving (linear) generalized
eigenproblems.

Input: Positive definite matrix B, factorized form of $\altmathcalb{L}(s)$,
size of Krylov subspace $m$, and initial vector $\textbf{v}_{1}$ .
Output: Arnoldi bases $\textbf{v}_{j}$ $(j=1,2,...,m+1)$ , upper Hessenberg
matrix $\textbf{H}_{m}$, and the residual factor $\beta$ .

1:Normalize:
$\textbf{v}_{1}\leftarrow\textbf{v}_{1}/\left\lVert\textbf{v}_{1}\right\rVert_{\\!\>\textbf{B}}$
, $\textbf{V}_{1}=\textbf{v}_{1}$

2:for all $j=1,2,...,m$, do

3:
$\textbf{r}=\textbf{S}\textbf{v}_{j}=\altmathcalb{L}(s)^{-1}\textbf{B}\textbf{v}_{j}$

4: $\textbf{h}=\textbf{V}_{j}^{T}\textbf{B}\textbf{r}$ ,
$h_{1:j,j}\leftarrow\textbf{h}$

5: $\textbf{z}=\textbf{r}-\textbf{V}_{j}\\!\>\textbf{h}$

6: $h_{j+1,j}=\left\lVert\\!\>\textbf{z}\\!\>\right\rVert_{\\!\>\textbf{B}}$

7: $\textbf{v}_{j+1}=\textbf{z}/h_{j+1,j}$ ,
$\textbf{V}_{j+1}\leftarrow\begin{bmatrix}\textbf{V}_{j}&\textbf{v}_{j+1}\end{bmatrix}$

8:end for

9:$\beta=h_{m+1,m}$

The Arnoldi vector $\textbf{v}_{m+1}$ is B-orthogonal with respect to the
columns of $\textbf{V}_{m}$ in the Gram–Schmidt sense, resulting in
${\textbf{V}_{m}^{T}\\!\>\textbf{B}\textbf{v}_{m+1}=0}$. Hence, the B-inner
product of $\textbf{V}_{m}$ premultiplied in (4.4) leads to the following
equation, noting that $\textbf{V}_{m}$ is B-orthogonal (i.e.,
${\textbf{V}_{m}^{T}\\!\>\textbf{B}\textbf{V}_{m}=\textbf{I}}$ ) :

$\displaystyle\textbf{V}_{m}^{T}\\!\>\textbf{B}\\!\>\textbf{S}\textbf{V}_{m}=\textbf{H}_{m}\,.$
(4.6)

In (4.6), $\textbf{H}_{m}$ is the B-orthogonal projection of S onto Krylov
subspace. Therefore, the $m$ eigenpairs of $\textbf{H}_{m}$ (referred to as
Ritz pairs) are the Rayleigh–Ritz approximation of the $m$ eigenpairs
${(\theta_{j},\tilde{\textbf{u}}_{j})}$ of the shift-and-invert problem (4.2).
Since $\textbf{H}_{m}$ is not symmetric, some eigenvalues of (4.2) may be
complex.

###### Remark 3

If A is Hermitian and ${s\in\mathbb{R}}$, the calculation of h in line 4 of
Algorithm 1 includes only the last basis vector (i.e.,
${h_{j,j}=\textbf{v}_{j}^{T}\textbf{B}\textbf{r}}$) and the Gram–Schmidt
orthogonalization (line 5) is only performed against the last two basis
vectors (i.e.,
${\textbf{z}=\textbf{r}-h_{j-1,j}\textbf{v}_{j-1}-h_{j,j}\textbf{v}_{j}}$). In
this case, $\textbf{H}_{m}$ converts to a symmetric tridiagonal matrix
resulting in real eigenvalues. Equation (4.4) is then referred to as the
Lanczos decomposition.

####  4.1.3 Restarting

When seeking a large number of eigenvalues ($N_{\rm eig}$), the norm of the
residual term in (4.4) grows, and the approximation quality of Ritz pairs
deteriorates. Thus, Krylov eigensolvers incorporate iterative restarting
mechanisms to maintain the residual norm lower than the desired tolerance
$\epsilon$. By setting $m$ slightly larger than $N_{\rm eig}$ , the
eigensolver keeps a set of good eigenpairs after each Arnoldi process — those
with residual norms below $\epsilon$. Then, a new $m$-step Arnoldi process
restarts in the next iteration, benefiting from the previously obtained
spectral approximation. The restarting process continues until we compute all
requested eigenpairs. Well-known restarting techniques are the thick-restart
Lanczos method [41], and its nonsymmetric equivalent for the Arnoldi case
referred to as the Krylov–Schur method [42, 43]. We use the latter in our
eigencomputations, but we omit details for the sake of brevity.

####  4.1.4 Spectrum slicing

The spectrum slicing technique is only applicable when all eigenvalues are
real [38, 44]. Thus, most quadratic eigensolvers cannot use this technique as
their spectrum may be complex333Under specific circumstances described in [25]
the spectrum slicing technique applies to QEPs. Still, we do not study these
cases here.. Nevertheless, we briefly review this algorithm to highlight a
significant difference between generalized and quadratic eigenproblems.

By definition, if all eigenvalues of an arbitrary matrix A are real, the
inertia of A is given by three integers ${\nu\\!\>(\textbf{A})}$,
$\zeta(\textbf{A})$, and ${\pi\\!\>(\textbf{A})}$ denoting the numbers of
negative, zero, and positive eigenvalues of A, respectively. For linear
Hermitian definite pencils, we consider an LDL Cholesky factorization of
$\altmathcalb{L}(s)$, which is
${\textbf{A}-s\\!\>\textbf{B}=\textbf{L}\textbf{D}\textbf{L}^{T}},$ being D a
diagonal matrix. Then, based on Sylvester’s law of inertia [45], one writes
${\nu\\!\>(\textbf{A}-s\\!\>\textbf{B})=\nu\\!\>(\textbf{D})}$, implying that
the number of eigenvalues of $\altmathcalb{L}(\lambda)$ smaller than $s$ is
equal to the number of negative eigenvalues of D. As a result, considering two
shifts $s_{k}$ and $s_{k+1}$, the spectrum subinterval $[s_{k},s_{k+1}]$ has
${\nu\\!\>(\textbf{A}-s_{k+1}\textbf{B})-\nu\\!\>(\textbf{A}-s_{k}\textbf{B})}$
eigenvalues including their multiplicities. The inertia information helps to
slice the spectrum efficiently into subintervals, leading to the efficient
spectrum slicing technique. When $N_{\rm eig}$ is large, spectrum slicing has
two main benefits. First, the eigensolution accuracy increases because we look
for eigenvalues in a smaller spectrum subinterval (i.e., eigenvalues are
closer to the shift). Second, spectrum slicing prevents the deterioration of
the convergence rate of the Arnoldi decomposition (4.4) because we project
onto a smaller Krylov subspace at each Arnoldi run. In this manner,
eigensolvers set the subspace size $m$ to be slightly higher than the number
of to-be-computed eigenvalues per shift (see, e.g., [6]). Each shift then
contains the corresponding Arnoldi recurrences and restarts described in
previous sections. Importantly, when dealing with the non-Hermitian or
indefinite pencils that result from the linearization of quadratic
eigenproblems, the lack of exact inertia (due to a complex spectrum) prevents
eigensolvers to perform this technique [46]. For such eigensystems, we select
only one target shift inside the spectrum region of interest.

###  4.2 Quadratic eigenproblems

We consider the quadratic eigenproblem,

$\displaystyle(\textbf{K}+\lambda\textbf{C}+\lambda^{2}\textbf{M})\\!\>\tilde{\textbf{u}}=\textbf{0}\,,$
(4.7)

characterized by the quadratic matrix pencil
${\altmathcalb{Q}(\lambda)\coloneqq\textbf{K}+\lambda\textbf{C}+\lambda^{2}\textbf{M}}$
. For ${N\times N}$ system matrices, this eigenproblem has $2N$ eigenpairs. In
many practical cases, $\altmathcalb{Q}(\lambda)$ is either Hermitian or
gyroscopic (c.f., our model problems in Section 2), which implies that the
eigenvalues may appear in a combination of complex conjugate pairs and real
values. Under special circumstances that all eigenvalues of a quadratic
Hermitian eigensystem may be purely real, we call $\altmathcalb{Q}(\lambda)$ a
hyperbolic pencil [47, 48]. But this particular case is not studied here.

####  4.2.1 Linearization

Eigensolvers perform different types of projection for computing the
eigenvalues of quadratic eigenproblems. Methods such as Jacobi–Davidson [49]
and second-order Arnoldi (SOAR) [50] directly project the quadratic
eigenproblem. However, a common way to solve (4.7) is to perform a
linearization by projecting an equivalent linearized eigenproblem of a doubled
size onto Krylov subspace (see, e.g., [20, 23, 27, 24]). Some linearizations
preserve the Hermitian structure of matrices [25]. Nevertheless, the obtained
linear pencil is indefinite in most cases. Herein, we use the following
linearization of (4.7), which applies to our two model problems:

$\displaystyle\Bigg{(}\underbrace{\begin{bmatrix}\textbf{0}&\textbf{I}\\\
-\textbf{K}&-\textbf{C}\end{bmatrix}}_{\textbf{A}}-\lambda\underbrace{\begin{bmatrix}\textbf{I}&\textbf{0}\\\
\textbf{0}&\textbf{M}\end{bmatrix}}_{\textbf{B}}\Bigg{)}\begin{bmatrix}\tilde{\textbf{u}}\\\
\lambda\tilde{\textbf{u}}\end{bmatrix}=\textbf{0}\,,$ (4.8)

where A and B are doubled-size matrices. After linearizing the quadratic
eigenproblems, we use the same numerical algorithms as the generalized
eigensystems. In the resulting linearized eigenproblem (4.8), matrix B is
positive definite, while matrix A is non-Hermitian. Thus, the Arnoldi
recurrence of Algorithm 1 holds. However, we deal with $2N\times 2N$ matrices
that will increase the computational cost of the eigensolution. We may resolve
this issue by exploiting the block structure of the linearization (4.8)
without explicitly constructing the doubled-size matrices A and B (see, e.g.,
[51, 25]). Considering

$\displaystyle\altmathcalb{L}(s)=\textbf{A}-s\\!\>\textbf{B}=\begin{bmatrix}-s\\!\>\textbf{I}&\textbf{I}\\\
-\textbf{K}&-\textbf{C}-s\\!\>\textbf{M}\end{bmatrix}\,,$ (4.9)

using the block factorization of $\altmathcalb{L}(s)$ and taking into account
the definition of the quadratic pencil $\altmathcalb{Q}(s)$, we write the
operator matrix as

$\displaystyle\textbf{S}=\altmathcalb{L}(s)^{-1}\textbf{B}=\begin{bmatrix}\altmathcalb{Q}(s)^{-1}&\textbf{0}\\\
\textbf{0}&\altmathcalb{Q}(s)^{-1}\end{bmatrix}\begin{bmatrix}-\textbf{C}-s\\!\>\textbf{M}&-\textbf{M}\\\
\textbf{K}&-s\\!\>\textbf{M}\end{bmatrix}\,.$ (4.10)

We then multiply the operator matrix with the Arnoldi basis vectors v (line 3
of Algorithm 1) by splitting v into two upper and lower $N\times 1$ blocks as
${\textbf{v}={\textstyle\footnotesize\begin{bmatrix}\textbf{v}_{u}\\\
\textbf{v}_{l}\end{bmatrix}}}$ . Hence, by writing
${\textbf{r}={\textstyle\footnotesize\begin{bmatrix}\textbf{r}_{u}\\\
\textbf{r}_{l}\end{bmatrix}}}$ in a similar way, we obtain

$\displaystyle\textbf{r}=\altmathcalb{L}(s)^{-1}\textbf{B}\textbf{v}=\begin{bmatrix}-\altmathcalb{Q}(s)^{-1}(s\\!\>\textbf{M}\textbf{v}_{u}+\textbf{C}\textbf{v}_{u}+\textbf{M}\textbf{v}_{l})\\\
s\\!\>\textbf{r}_{u}+\textbf{v}_{u}\end{bmatrix}\\!\>.$ (4.11)

In this context, we write every B-norm calculation as:

$\displaystyle\left\lVert\textbf{v}\right\rVert_{\\!\>\textbf{B}}^{2}=\textbf{v}^{T}\textbf{B}\textbf{v}=\begin{bmatrix}\textbf{v}_{u}^{T}&\textbf{v}_{l}^{T}\end{bmatrix}\begin{bmatrix}\textbf{I}&\textbf{0}\\\
\textbf{0}&\textbf{M}\end{bmatrix}\begin{bmatrix}\textbf{v}_{u}\\\
\textbf{v}_{l}\end{bmatrix}=\textbf{v}_{u}^{T}\textbf{v}_{u}+\textbf{v}_{l}^{T}\textbf{M}\textbf{v}_{l}\\!\>.$
(4.12)

Algorithm 2 Numerical eigenanalysis procedure for solving quadratic
eigenproblems.

Input: ${N\times N}$ system matrices K, C and M, target shift $s$, requested
number of eigenvalues $N_{\rm eig}$ , size of Krylov subspace ${m>N_{\rm
eig}}$ , initial vector $\textbf{v}_{1}$ , and residual tolerance
$\epsilon$.
Output: eigenpairs $\lambda$ and $\tilde{\textbf{u}}$ .

1:Normalize:
$\textbf{v}_{1}\leftarrow\textbf{v}_{1}/\left\lVert\textbf{v}_{1}\right\rVert_{\\!\>\textbf{B}}$
, $\textbf{V}_{1}=\textbf{v}_{1}$

2:Factorize $\altmathcalb{Q}(s)$

3:$\altmathcal{N}_{\rm Iter}=0$ /* a counter for total number of iterations */

4:${R=10^{\\!\>6}}$ /* setting a big value for residual norm to get into
iterations */

5:while $R>\epsilon$, do

6: $\altmathcal{N}_{\rm Iter}\leftarrow\altmathcal{N}_{\rm Iter}+1$

7: for all $j=1,2,...,m$, do /* $m$-step Arnoldi recurrence */

8:
$\textbf{r}={\textstyle\footnotesize\begin{bmatrix}-\altmathcalb{Q}(s)^{-1}\big{(}s\\!\>\textbf{M}\textbf{v}_{j,u}+\textbf{C}\textbf{v}_{j,u}+\textbf{M}\textbf{v}_{j,l}\big{)}\\\
s\\!\>\textbf{r}_{u}+\textbf{v}_{j,u}\end{bmatrix}}$

9:
$\textbf{h}=\textbf{V}_{j}^{T}\textbf{B}\textbf{r}={\textstyle\footnotesize\begin{bmatrix}\textbf{v}_{1,u}^{T}\textbf{r}_{u}+\textbf{v}_{1,l}^{T}\\!\>\textbf{M}\textbf{r}_{l}\\\\[-3.0pt]
\vdots\\\
\textbf{v}_{j,u}^{T}\textbf{r}_{u}+\textbf{v}_{j,l}^{T}\\!\>\textbf{M}\textbf{r}_{l}\end{bmatrix}}$
, $h_{1:j,j}\leftarrow\textbf{h}$

10: $\textbf{z}=\textbf{r}-\textbf{V}_{j}\\!\>\textbf{h}$

11:
$h_{j+1,j}=\left\lVert\\!\>\textbf{z}\\!\>\right\rVert_{\\!\>\textbf{B}}=\big{(}\textbf{z}_{u}^{T}\textbf{z}_{u}+\textbf{z}_{l}^{T}\textbf{M}\textbf{z}_{l}\big{)}^{1/2}$

12: $\textbf{v}_{j+1}=\textbf{z}/h_{j+1,j}$ ,
$\textbf{V}_{j+1}\leftarrow\begin{bmatrix}\textbf{V}_{j}&\textbf{v}_{j+1}\end{bmatrix}$

13: end for

14: $\beta=h_{m+1,m}$

15: Construct matrix $\textbf{H}_{m}$ and compute Ritz pairs

16: Compute numerical eigenpairs $\lambda_{j}$ and $\tilde{\textbf{u}}_{j}$
$(j=1,2,...,m)$

17: Compute largest residual norm $R$

18: Keep adequate numbers of eigenpairs and compute initial values for restart

19:end while

This block representation allows us to perform factorization and matrix–vector
operations of the Arnoldi recurrence using ${N\times N}$ system matrices.
Other steps of solving quadratic eigenproblems exploit the same block-product
structure. Algorithm 2 presents a pseudo-code describing the numerical
eigenanalysis procedure for solving quadratic eigenproblems. However, the
entire process might be more complicated. For instance, the size of Krylov
subspace and the adequate number of Arnoldi bases to keep before restart play
essential roles in the efficiency of the eigensolution. These variables mainly
depend on the size of the eigensystem, the total number of requested
eigenpairs, and the selected eigensolver package.

##  5 Eigencomputation cost

###  5.1 General overview

The total eigencomputation cost of quadratic eigenproblems is mostly governed
by the following operation sets

  * 1.

LU factorization (Fact) of ${\altmathcalb{Q}(s)}$;

  * 2.

Forward/backward (FB) eliminations (i.e., multiplications of LU factors by
vectors);

  * 3.

Matrix–vector (MV) and vector–vector (VV) products in the sense of Krylov
projection (see Algorithm 2).

We consider other costs (e.g., those related to creating system matrices) as
lower-order terms and exclude them from our cost estimation. Thus, we write

$\displaystyle\altmathcal{C}\approx\altmathcal{N}_{\rm
Fact}\\!\>\altmathcal{C}_{\rm Fact}+\altmathcal{N}_{\rm
FB}\\!\>\altmathcal{C}_{\rm FB}+\altmathcal{N}_{\rm
MV}\\!\>\altmathcal{C}_{\rm MV}+\altmathcal{N}_{\rm
VV}\\!\>\altmathcal{C}_{\rm VV}\,,$ (5.1)

where with $\altmathcal{C}$ and $\altmathcal{N}$, we refer to the cost and
number of operations, respectively.

For linear Hermitian definite pencils, the spectrum slicing entails performing
multiple factorizations (one per shift). At the same time, thanks to a smaller
Krylov subspace, fewer matrix–vector operations are required (see Section
4.1.4). Additionally, in such eigensystems, only the last two basis vectors
participate in each recurrence step (see Remark 3), implying that the cost of
vector–vector products is negligible. As a result, for moderate to large
generalized eigensystems (characterized by linear Hermitian definite pencils),
matrix factorization governs the total eigencomputation cost (see [6]). On the
other hand, for quadratic eigenproblems with non-Hermitian (or indefinite)
linearized pencils, we perform one shift and, therefore, one factorization for
all requested eigenpairs, resulting in a larger Krylov subspace (the subspace
size grows as $N_{\rm eig}$ increases). Moreover, a full basis matrix V
participates in the Arnoldi recurrence. Thus, when solving moderately sized
problems, the other three operations in (5.1) govern the eigencomputation
cost. Their contribution to the total cost is remarkable when we look for a
large number of eigenpairs. Nevertheless, the matrix factorization is
asymptotically the most expensive operation as the problem size grows.

###  5.2 Eigencomputation cost per operation

####  5.2.1 Matrix factorization

As we state above, solving quadratic eigenproblems entails one LU
factorization of the quadratic pencil ${\altmathcalb{Q}(s)}$ for the target
shift $s$ inside the eigenspectrum region of interest (i.e.,
${\altmathcal{N}_{\rm Fact}=1}$, see line 2 of Algorithm 2). To determine the
factorization cost, we use the theoretical estimates of IGA and rIGA
discretizations for H(curl) and H(div) spaces [3]. Thus, for 2D problems, we
write

$\displaystyle\altmathcal{C}_{{\rm Fact}\\!\>,\\!\>{\rm IGA}}$
$\displaystyle\approx\altmathcal{O}\big{(}N^{1.5}(p-1)^{3}\big{)}\,,$ (5.2)
$\displaystyle\altmathcal{C}_{{\rm Fact}\\!\>,\\!\>{\rm rIGA}}$
$\displaystyle\approx\altmathcal{O}\big{(}N^{1.5}(p-1)\big{)}\,.$ (5.3)

The number of degrees of freedom of an rIGA system is slightly higher than its
IGA counterpart (a.k.a., ${N_{\rm rIGA}\approx N_{\rm IGA}}$). Thus, we
consider ${N=N_{\rm IGA}}$ in both (5.2) and (5.3) and obtain a cost
improvement of ${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$ for large problem
sizes when using rIGA discretization with optimal macroelements containing 16
elements in each direction. For smaller eigensystems, computational savings
above ${\altmathcal{O}(p-1)}$ are hardly obtainable (see [3]).

####  5.2.2 Forward/backward elimination

The number of forward/backward eliminations is equal to the size of Krylov
subspace, $m$, multiplied by the number of iterations, $\altmathcal{N}_{\rm
Iter}$ (i.e., ${\altmathcal{N}_{\rm FB}=m\\!\>\altmathcal{N}_{\rm Iter}}$ ,
see line 8 of Algorithm 2 in computing the upper block of r). The
forward/backward elimination cost is proportional to the number of nonzero
terms of the LU factors, which we estimate for 2D IGA and rIGA discretizations
as follows (see [1, 6]) :

$\displaystyle\altmathcal{C}_{{\rm FB}\\!\>,\\!\>{\rm IGA}}$
$\displaystyle\approx\altmathcal{O}\big{(}N(p-1)^{2}\big{)}\,,$ (5.4)
$\displaystyle\altmathcal{C}_{{\rm FB}\\!\>,\\!\>{\rm rIGA}}$
$\displaystyle\approx\altmathcal{O}\big{(}N(p-1)\big{)}\,.$ (5.5)

Again, the cost reduction of ${\altmathcal{O}(p-1)}$ with respect to IGA in
(5.5) is asymptotically apparent only for large problem sizes and when
employing macroelements with 16 elements in each direction.

####  5.2.3 Matrix–vector multiplication

The quadratic eigenproblem solver calls the matrix–vector operator in lines 8,
9 and 11 of Algorithm 2. Thus, one obtains ${\altmathcal{N}_{\rm
MV}=5m\\!\>\altmathcal{N}_{\rm Iter}}$ for the participation of either M or C
in the process. This numerical operation has a cost proportional to the number
of nonzero entries, $N_{nz}$ , of system matrices — that is related to the
sum of interactions of each basis function with all other bases [52]. For 2D
eigensystems, we write

$\displaystyle\altmathcal{C}_{{\rm MV}\\!\>,\\!\>{\rm IGA}}$
$\displaystyle=\altmathcal{O}\big{(}N_{nz}(\textbf{M}_{\rm
IGA})\big{)}\approx\altmathcal{O}(N_{\rm IGA}p^{2})\,,$ (5.6)
$\displaystyle\altmathcal{C}_{{\rm MV}\\!\>,\\!\>{\rm rIGA}}$
$\displaystyle=\altmathcal{O}\big{(}N_{nz}(\textbf{M}_{\rm
rIGA})\big{)}\approx\altmathcal{O}(N_{\rm rIGA}p^{2})\,.$ (5.7)

An rIGA-discretized system has a slightly higher number of degrees of freedom.
Thus, matrix–vector products suffer a slight degradation when using rIGA. We
distinguish between $N_{\rm IGA}$ and $N_{\rm rIGA}$ in (5.6) and (5.7) since
it is the only difference in the cost of this numerical operation with respect
to the employed discretization.

####  5.2.4 Vector–vector multiplication

The vector–vector product is theoretically the cheapest numerical operation of
eigenanalysis. It has a computational complexity proportional to the number of
entries of vectors, $N$, which is significantly smaller than
$N_{nz}(\textbf{M})$. However, when solving quadratic eigenproblems, we
perform a large number of vector–vector products in the sense of Gram–Schmidt
orthogonalization (lines 8 –11 of Algorithm 2). In particular, at the
${j\\!\>}$-th Arnoldi step ${(j=1,2,...,m)}$, the participation of
$\textbf{V}_{j}$ incorporates with $2j$ multiplications of vectors of size $N$
in line 9, and $2N$ products of vectors of size $j$ in line 10. At each
iteration, the eigensolver repeats these lines $m$ times during the Arnoldi
process. Thus, one obtains ${\altmathcal{N}_{\rm VV}\approx
2m^{2}\altmathcal{N}_{\rm Iter}}$ as the total number of multiplications of
vectors of size $N$. As a result, the vector–vector products have a comparable
cost as matrix–vector operations when we seek to compute a large number of
eigenvalues (i.e., a large $m$). In contrast to other operations, the cost of
each vector–vector product is not directly related to the polynomial degree
$p$ as it only depends on $N$ :

$\displaystyle\altmathcal{C}_{{\rm VV}\\!\>,\\!\>{\rm IGA}}$
$\displaystyle=\altmathcal{O}(N_{\rm IGA})\,,$ (5.8)
$\displaystyle\altmathcal{C}_{{\rm VV}\\!\>,\\!\>{\rm rIGA}}$
$\displaystyle=\altmathcal{O}(N_{\rm rIGA})\,.$ (5.9)

###  5.3 Summary of eigencomputation cost

Equations (5.2) and (5.3) show that matrix factorization is the most expensive
operation when $N$ is large. The other three operations in (5.1) are not
inherently as expensive as the matrix factorization. However, the eigensolver
calls them many times during the Krylov projection. Thus, they become decisive
when the number of requested eigenvalues $N_{\rm eig}$ is large. Since the
subspace size $m$ is slightly higher than $N_{\rm eig}$ , we assume
${m\approx N_{\rm eig}}$ and summarize the contribution of each numerical
operation to the total eigencomputation cost as tabulated in Table 1.

Table 1: Most expensive operations of the quadratic eigensolver and their contribution to the total eigencomputation cost for 2D eigensystems. Cost improvements of rIGA are under the assumption of being $N$ sufficiently large and using optimal macroelements with 16 elements in each direction. Unless stated otherwise, ${N=N_{\rm IGA}}$ . Numerical operation |  | LU factorization | FB elimination | Mat–vec product | Vec–vec product
---|---|---|---|---|---
Number of times the operation is called |  | 1 | $N_{\rm eig}\\!\>\altmathcal{N}_{\rm Iter}$ | $5N_{\rm eig}\\!\>\altmathcal{N}_{\rm Iter}$ | $2N_{\rm eig}^{2}\\!\>\altmathcal{N}_{\rm Iter}$
Cost of performing one operation | IGA | $\altmathcal{O}(N^{1.5}p^{3})$ | $\altmathcal{O}(Np^{2})$ | $\altmathcal{O}(N_{\rm IGA}p^{2})$ | $\altmathcal{O}(N_{\rm IGA})$
| rIGA | $\altmathcal{O}(N^{1.5}p)$ | $\altmathcal{O}(Np)$ | $\altmathcal{O}(N_{\rm rIGA}p^{2})$ | $\altmathcal{O}(N_{\rm rIGA})$
Total cost of performing the operation | IGA | $\altmathcal{O}(N^{1.5}p^{3})$ | $\altmathcal{O}(N_{\rm eig}Np^{2})$ | $\altmathcal{O}(N_{\rm eig}N_{\rm IGA}p^{2})$ | $\altmathcal{O}(N_{\rm eig}^{2}N_{\rm IGA})$
| rIGA | $\altmathcal{O}(N^{1.5}p)$ | $\altmathcal{O}(N_{\rm eig}Np)$ | $\altmathcal{O}(N_{\rm eig}N_{\rm rIGA}p^{2})$ | $\altmathcal{O}(N_{\rm eig}^{2}N_{\rm rIGA})$
Improvement/degradation of |  | Improved by | Improved by | Degraded by | Degraded by
performing the operation in rIGA |  | $\altmathcal{O}(p^{2})$ | $\altmathcal{O}(p)$ | $N_{\rm rIGA}/N_{\rm IGA}$ | $N_{\rm rIGA}/N_{\rm IGA}$

##  6 Numerical results

We report numerical results of our eigencomputations when using maximum-
continuity and refined isogeometric analyses to solve the electromagnetic and
vibroacoustic eigenproblems described in Section 2. We study the numerical
efficiency and accuracy versus cost for different IGA and rIGA
discretizations. Before proceeding with results, we provide implementation
details in the following subsection.

###  6.1 Implementation details

We discretize our model problems using PetIGA-MF [53], which is a multifield
extension of PetIGA [54] — a high-performance isogeometric analysis framework
based on PETSc [55]. PetIGA-MF, utilized in many scientific and engineering
applications (see, e.g., [56, 57, 58, 59, 60, 3, 61, 7]), uses different
spaces for each field of interest and employs data management libraries to
condense the data of multiple fields in a single object, thus simplifying the
discretization construction.

We also use SLEPc, the scalable library for eigenvalue problem computations
[62], for performing the eigenanalysis. SLEPc, implemented to solve different
eigenproblems types (see, e.g., [44, 63, 24, 64, 65, 66]), allows us to employ
the quadratic eigenproblem solver incorporating the shift-and-invert spectral
transformation, linearization, the Arnoldi decomposition, and the Krylov–Schur
restarting methods.

We use multifrontal direct solver MUMPS [67] to perform matrix factorization
and the forward/backward eliminations. We employ the sequential version of
MUMPS, which runs on a single thread (core). We also use the automatic matrix
reordering algorithm provided by METIS [35].

###  6.2 Case studies

####  6.2.1 Electromagnetic wave propagation in a three-layer Earth model

Figure 4: A three-layer Earth model as the physical domain of the
electromagnetic eigenproblem. Each layer has electrical conductivities (given
in ${\rm S/m}$) along horizontal and vertical directions.

We investigate the electromagnetic wave propagation eigenproblem. In
particular, we consider a three-layer heterogeneous Earth model as a common
case study in geosteering applications (see, e.g., [7, 68]). This 2D Earth
model is essential in studying the full 3D wave propagation eigenproblem in
transversely isotropic media [69]. As Fig. 4 shows, in each layer, we consider
electrical conductivities along horizontal and vertical directions given by
the conductivity matrix $\boldsymbol{\sigma}$ as follows:

$\displaystyle\boldsymbol{\sigma}\coloneqq\begin{bmatrix}\sigma_{x}&0\\\
0&\sigma_{y}\end{bmatrix}\\!\>.$ (6.1)

Fig. 5 shows a few approximate eigenvalues of the quadratic eigenproblem
(2.13) arising in the electromagnetic wave propagation through our three-layer
Earth model. For the sake of simplicity and to have a better visual
representation of eigenvalues, we assume ${\mu=1}$ , and ${\varepsilon=1}$ .
We represent the approximate eigenfunctions associated with two arbitrary
eigenmodes in Fig. 6 to see how the electromagnetic wave looks like in
different layers. The electrical conductivity (i.e., the inverse of
resistivity) performs as a proportional damping. Thus, eigenfunctions show a
non-decaying oscillatory pattern within each layer. The decaying effect of
material conductivities, however, appears in the time response that is not
studied here (see [22] for further details).

Figure 5: Approximate eigenvalues for electromagnetic wave propagation
eigenproblem of the three-layer Earth model; $\omega_{h}$ values are given in
${\rm rad/s}$.

(a) Approximate eigenfunctions associated with
$\omega_{h}=7.024+0.282\\!\>{\rm i\\!\>}~{}{\rm rad/s}$

(b) Approximate eigenfunctions associated with
$\omega_{h}=19.200+0.193\\!\>{\rm i\\!\>}~{}{\rm rad/s}$

Figure 6: Two sample approximate eigenfunctions of the electromagnetic
eigenproblem of Fig. 4; $\textbf{E}_{h}$ values are given in ${\rm V/m}$.

####  6.2.2 Acoustic fluid in a cavity with one absorbing wall

A practical application of the acoustic fluid in a cavity with absorbing walls
occurs in the problem of decreasing noise level in, for example, a vehicle
cabin. Absorbing walls, commonly covered by a thin layer of a viscoelastic
material, can dissipate the fluid’s acoustic energy. In here, we consider a 2D
rectangular cavity with one absorbing wall, as Fig. 7 depicts. We assume the
cavity is filled with air ${(\,\rho=1~{}{\rm kg/m^{3}})}$ , the acoustic
speed is ${c=340~{}{\rm m/s}}$ , and the viscoelastic material has the
impedance parameters ${\alpha=5\times 10^{4}~{}{\rm N/m^{3}}}$ and
${\beta=200~{}{\rm N\cdot s/m^{3}}}$.

Figure 7: Physical domain of the vibroacoustic eigenproblem. We assume the
cavity is filled with air and has one absorbing wall.

Fig. 8 shows a few approximate eigenvalues of the second case study. Again,
the eigenvalues form a combination of conjugate complex pairs and purely real
values. The negative real part of the eigenvalues connotes the decaying
behavior of vibration modes of the fluid. Purely real eigenvalues,
corresponding to overdamped modes, theoretically have two accumulation points
at $-\frac{\alpha}{\beta}$ and $-\frac{2\alpha}{\beta}$ (see [70, 23] for more
details). More precisely, the eigenvalues exist in two separate branches,
namely ${-\frac{2\alpha}{\beta}<\lambda_{h}<-\frac{\alpha}{\beta}}$ and
${\lambda_{h}<-\frac{2\alpha}{\beta}}$ , which correspond to two possible
solutions of the dispersion equations (see Section 6.4). Fig. 10 illustrates
two eigenfunctions of two distinct eigenmodes in two different branches. The
dissipation effect of the absorbing boundary causes the solution to decay in
the vertical direction while it is oscillatory along the horizontal one.

Figure 8: Approximate eigenvalues for quadratic eigenproblem arising from the
acoustic fluid motion in a rectangular cavity with one absorbing wall; we
represent only the first accumulation point of purely real eigenvalues at
${-\frac{\alpha}{\beta}=-250~{}{\rm rad/s}}$.

(a) Numerical eigenfunctions associated with $\lambda_{h}=-341.80~{}{\rm
rad/s}$

(b) Numerical eigenfunctions associated with $\lambda_{h}=-3480.48~{}{\rm
rad/s}$

Figure 10: Two arbitrary normalized eigenfunctions corresponding to two
distinct purely real eigenvalues in two different branches of the
vibroacoustic case study.

###### Remark 4

When the stiffness, mass, and damping matrices have widely different norms,
which could be the case in the current vibroacoustic eigenproblem, the
solution of the linearized eigensystem (4.8) may suffer a significant
numerical error (c.f., [71]). We overcome this issue, as recommended in [71],
by solving a scaled eigenproblem:

$\displaystyle\big{(}\textbf{K}+\gamma\\!\>\textbf{C}_{\varsigma}+\gamma^{\\!\>2}\textbf{M}_{\varsigma}\big{)}\\!\>\tilde{\textbf{u}}=\textbf{0}\,,$
(6.2)

in place of the original one, where ${\gamma=\lambda_{h}/\varsigma}$ ,
${\textbf{C}_{\varsigma}=\varsigma\\!\>\textbf{C}}$ , and
${\textbf{M}_{\varsigma}=\varsigma^{\\!\>2}\textbf{M}}$ . In this setup, the
scaling factor $\varsigma$ is

$\displaystyle\varsigma\coloneqq\sqrt{\dfrac{\left\lVert\textbf{K}\right\rVert_{\infty}}{\left\lVert\textbf{M}\right\rVert_{\infty}}}\,,$
(6.3)

where the infinity norm is obtained as the maximum absolute value of
respective matrices.

###  6.3 Computational efficiency tests

We report the computational cost of the eigensolution of our two cases
studies. Since the computational cost of solving quadratic eigensystems
discretized by either H(curl) or H(div) spaces are approximately the same, we
only present the costs of solving the eigenproblem of electromagnetic wave
propagation. We test different mesh sizes with ${n_{e}=2^{s}}$
${(s=5,6,...,11)}$ elements in each direction, and different polynomial
degrees of B-spline bases, namely ${p=3,4,5,6}$. We assess the computational
efficiency of the rIGA framework by considering an optimal macroelement size
for rIGA discretization; that is, each macroelement consists of a ${16\times
16}$ grid (c.f., [1, 3, 6]).

Fig. 11 describes the contribution of each of the most expensive operations to
the total eigencomputation cost when using maximum-continuity and refined IGA
frameworks to approximate ${N_{\rm eig}=500}$ eigenpairs. We also investigate
the cost improvement (or degradation) rates of each operation under the
employment of rIGA by plotting the relative cost ${\altmathcal{C}_{\rm
IGA}/\altmathcal{C}_{\rm rIGA}}$ against the mesh size in Fig. 12. For large
problems, results indicate an improvement in the factorization cost close to
${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$ and of approximately
${\altmathcal{O}(p-1)}$ for forward/backward eliminations. There is a slight
degradation in the cost of matrix–vector and vector–vector multiplications due
to the slight increase of the system size in rIGA. In summary, the total
observed cost saving for the entire eigensolution is up to
${\altmathcal{O}(p-1)}$ for large domains (see Table 2). If the problem size
is large, the total computational cost is governed only by matrix
factorization. Therefore, we predict the total time improvements of up to
${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$. To observe this scaling, we would
need larger computational resources. On the other hand, for small problems
(e.g., systems with ${n_{e}\leq 128}$), the overall cost of maximum-continuity
IGA is comparable (or even smaller) to that of rIGA, which occurs when the
matrix factorization cost is a small fraction of the total cost. Finally, the
number of requested eigenpairs affects our cost improvements (see Section
5.3). However, in the case of a very large $N_{\rm eig}$ , the eigenvalues
may be located far from the target shift, thus, decreasing the accuracy of the
spectral approximation. Consequently, a reasonable number of eigenpairs should
be sought.

(a) $p=3$

(b) $p=4$

(c) $p=5$

(d) $p=6$

Figure 11: Computational cost (FLOPs number) of different numerical operations
employed by the eigensolver when using IGA and rIGA discretizations. We report
the cost of finding 500 eigenpairs of quadratic eigenproblems discretized by
H(curl) and H(div) spaces. Due to highly demanding computations, for
${n_{e}=2048}$, we only test up to ${p=5}$.

Figure 12: Cost improvement (or degradation) rates of different numerical operations of the eigenanalysis for rIGA when finding 500 eigenpairs of quadratic eigenproblems discretized by H(curl) and H(div) spaces.  Table 2: Improvement/degradation rates of different numerical operations of the eigenanalysis when using rIGA discretizations. We report the results associated with ${n_{e}=1024,2048}$ in Fig. 12 to show the asymptotic improvement factor of rIGA. |  | Factorization | FB elimination | Mat–vec | Vec–vec | Total cost
---|---|---|---|---|---|---
Degree |  $n_{e}$ | improved by | improved by | degraded by | degraded by | improved by
$p=3$ |  1024 | 3.634 | 1.937 | 0.923 | 0.940 | 1.376
2048 | 3.696 | 1.986 | 0.921 | 0.940 | 1.602
$p=4$ |  1024 | 7.721 | 2.825 | 0.867 | 0.886 | 1.839
2048 | 8.145 | 2.957 | 0.865 | 0.886 | 2.454
$p=5$ |  1024 | 11.634 | 3.622 | 0.811 | 0.836 | 2.303
2048 | 14.775 | 4.124 | 0.810 | 0.836 | 3.549
$p=6$ |  1024 | 18.203 | 4.097 | 0.765 | 0.791 | 2.811
2048 | *** | *** | *** | *** | ***

###  6.4 Accuracy versus cost

We present the results of accuracy-versus-cost tests when using different IGA
and rIGA discretizations to solve eigenproblems arising in electromagnetics
and vibroacoustics. In [18, 6], the authors show for linear Hermitian definite
pencils that, other than for the outliers, both IGA and rIGA frameworks result
in almost the same accuracy in the eigensolution. In the following, we
consider one arbitrary eigenmode for each case study and compare the
eigenvalue and eigenfunction errors. We show that rIGA delivers almost the
same accuracy than IGA, but with a lower computational effort.

####  6.4.1 Electromagnetic wave propagation eigenproblem

The analytic eigensolution of (2.3) in conductive media is unavailable. Thus,
we consider ${\boldsymbol{\sigma}=\textbf{0}}$ and study the accuracy of
approximate eigensolution. The analytic eigenvalues and eigenfunctions of the
electromagnetic wave in a non-conductive unit square are expressed in a tensor
form as follows (see, e.g., [10]):

$\displaystyle\lambda_{ij}$
$\displaystyle=\pi^{\\!\>2}(i^{\\!\>2}+j^{\\!\>2})\,,$ (6.4)
$\displaystyle\textbf{E}_{ij}$
$\displaystyle=\dfrac{2}{\sqrt{i^{\\!\>2}+j^{\\!\>2}}}\begin{bmatrix}-j\cos(i\pi
x)\sin(j\pi y)\\\\[2.0pt] ~{}i\sin(i\pi x)\cos(j\pi y)\end{bmatrix}\\!\>.$
(6.5)

Figs. 13(a) and 13(b) show the total computational cost per eigenmode versus
the eigenvalue error ${(\lambda_{h}-\lambda)/\lambda}$ and eigenfunction
$L^{2}$-norm error
${\left\lVert\\!\>\textbf{E}_{h}-\textbf{E}\\!\>\right\rVert_{\\!\>L^{2}}}$ ,
respectively. We consider the analytic eigenvalue
${\lambda_{23,23}=10442.04~{}{\rm rad/s}}$ and use IGA and rIGA
discretizations with different polynomial degrees and mesh sizes. Both figures
confirm that under the employment of rIGA, we obtain almost the same
approximation error as the maximum-continuity IGA with a lower computational
effort. Here, in the absence of conductivity, the quadratic eigenproblem is
converted to a definite generalized one. Thus, we follow the notations
described in Remark 2 and use $\lambda$ in reference to eigenvalues of the
electromagnetic wave. Additionally, we use the generalized eigensolver to
solve the eigenproblem, thus, obtaining higher improvements in the number of
FLOPs when using rIGA.

(a) Eigenvalue error: $(\lambda_{h}-\lambda)/\lambda$

(b) Eigenfunction $L^{2}$-norm error:
$\left\lVert\\!\>\textbf{E}_{h}-\textbf{E}\\!\>\right\rVert_{\\!\>L^{2}}$

Figure 13: Accuracy versus eigencomputation cost for the electromagnetic
eigenvalue problem in a non-conductive unit square. We compare the results of
different IGA and rIGA discretizations when approximating the eigenpair
${(\lambda_{23,23},\textbf{E}_{23,23})}$.

####  6.4.2 Acoustic fluid eigenproblem

We obtain the analytical vibration modes of the acoustic fluid inside a 2D
rectangular cavity with one absorbing wall (Fig. 7) by separation of variables
as follows (see [70]):

$\displaystyle\textbf{u}_{j}=\dfrac{-1}{\rho\lambda_{j}^{2}}\nabla\left(\cos\Big{(}\dfrac{j\pi
x}{a}\Big{)}\cosh(\eta_{j}y)\right)=\dfrac{-1}{\rho\lambda_{j}^{2}}\begin{bmatrix}-\dfrac{j\pi}{a}\sin\Big{(}\dfrac{j\pi
x}{a}\Big{)}\cosh(\eta_{j}y)\\\\[8.0pt] ~{}\eta_{j}\cos\Big{(}\dfrac{j\pi
x}{a}\Big{)}\sinh(\eta_{j}y)\end{bmatrix}\\!\>,$ (6.6)

where eigenvalues $\lambda_{j}$ and auxiliary variables
${\eta_{j}\in\mathbb{C}}$ are solutions of the dispersion equations:

$\displaystyle\eta_{j}^{\\!\>2}$
$\displaystyle=\dfrac{\lambda_{j}^{2}}{c^{\\!\>2}}+\dfrac{\pi^{\\!\>2}}{a^{\\!\>2}}j^{\\!\>2}\,,$
(6.7) $\displaystyle\eta_{j}\tanh(\eta_{j}\\!\>b)$
$\displaystyle=-\dfrac{\rho\lambda_{j}^{2}}{\alpha+\lambda_{j}\\!\>\beta}\,.$
(6.8)

Herein, we seek purely real eigenvalues. As Section 6.2.2 states, for each
$j$, the eigenvalues have two branches that accumulate at
$-\frac{\alpha}{\beta}$ and $-\frac{2\alpha}{\beta}$. Figs. 14(a) and 14(b)
show the eigenvalue error ${(\lambda_{h}-\lambda)/\lambda}$ and eigenfunction
$L^{2}$-norm error
${\left\lVert\\!\>\textbf{u}_{h}-\textbf{u}\\!\>\right\rVert_{\\!\>L^{2}}}$
versus the total eigencomputation cost per eigenmode when using different IGA
and rIGA discretizations. The results are associated with
${\lambda_{10}=-7377.39~{}{\rm rad/s}}$ on the second branch (i.e.,
${\lambda<-\frac{2\alpha}{\beta}}$). Again, the plots confirm that the
implementation of rIGA reduces the computational cost for almost the same
numerical error.

(a) Eigenvalue error: $(\lambda_{h}-\lambda)/\lambda$

(b) Eigenfunction $L^{2}$-norm error:
$\left\lVert\\!\>\textbf{u}_{h}-\textbf{u}\\!\>\right\rVert_{\\!\>L^{2}}$

Figure 14: Accuracy versus eigencomputation cost for the quadratic
eigenproblem arising in vibroacoustics. We compare the results of different
IGA and rIGA discretizations when approximating the eigenpair
${(\lambda_{10},\textbf{u}_{10})}$ on the second eigenbranch of the
rectangular cavity of Fig. 7.

##  7 Conclusions

We use refined isogeometric analysis (rIGA) discretizations to solve quadratic
eigenvalue problems (QEPs). We compare the computational performance of rIGA
versus that of maximum-continuity IGA when solving quadratic eigensystems
arising in 2D multifield problems discretized by H(curl) and H(div) spaces.
For large problem sizes, the most expensive numerical operation of the
eigenanalysis is matrix factorization, followed by matrix–vector and
vector–vector operations in the sense of Krylov projections. We estimate that
the LU factorization is asymptotically
${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$ times faster with an rIGA
discretization compared to IGA. As a result, we can theoretically reach an
improvement of ${\altmathcal{O}\big{(}(p-1)^{2}\big{)}}$ in the total
eigencomputation cost of vector-valued multifield eigenproblems — while it is
expected to be ${\altmathcal{O}(p^{2})}$ in scalar-valued problems in $H^{1}$.
Our numerical tests show that the computational savings associated with the
implementation of rIGA tend to ${\altmathcal{O}(p-1)}$, as we need to analyze
larger systems to observe the asymptotic behavior. For small problems, the
improvement hardly reaches the expected rates. Thus, we suggest using the
maximum-continuity IGA discretization only for small problems. Finally,
quadratic eigensystems have complex eigenvalues in most occasions, forcing the
eigensolver to project the problem onto a larger Krylov subspace (as compared
to the definite generalized eigenproblems). Thus, more matrix–vector and
vector–vector operations are required, resulting in a deterioration of the
improvement rate of rIGA versus IGA when solving quadratic eigenproblems.

## Acknowledgment

This work has received funding from the European Union’s Horizon 2020 research
and innovation program under the Marie Sklodowska-Curie grant agreement No.
777778 (MATHROCKS); the European Regional Development Fund (ERDF) through the
Interreg V-A Spain-France-Andorra program POCTEFA 2014-2020 Project PIXIL
(EFA362/19); the Spanish Ministry of Science and Innovation projects with
references PID2019-108111RB-I00 (FEDER/AEI) and PDC2021-121093-I00, the “BCAM
Severo Ochoa” accreditation of excellence (SEV-2017-0718); and the Basque
Government through the BERC 2018-2021 program, the three Elkartek projects
3KIA (KK-2020/00049), EXPERTIA (KK-2021/00048), and SIGZE (KK-2021/00095), and
the Consolidated Research Group MATHMODE (IT1294-19) given by the Department
of Education. This publication was also made possible in part by the
Professorial Chair in Computational Geoscience at Curtin University and the
Deep Earth Imaging Enterprise Future Science Platforms of the Commonwealth
Scientific Industrial Research Organisation, CSIRO, of Australia. The Curtin
Corrosion Centre and the Curtin Institute for Computation kindly provide
ongoing support. The authors also acknowledge the computer resources and
technical support provided by Barcelona Supercomputing Center through the
MareNostrum4 cluster (activity IDs IM-2020-3-0009 and IM-2021-2-0015).

## References

## References

  * Garcia et al. [2017] D. Garcia, D. Pardo, L. Dalcin, M. Paszyński, N. Collier, V. M. Calo, The value of continuity: Refined isogeometric analysis and fast direct solvers, Computer Methods in Applied Mechanics and Engineering 316 (2017) 586–605, doi:10.1016/j.cma.2016.08.017.
  * Hughes et al. [2005] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (39-41) (2005) 4135–4195, doi:10.1016/j.cma.2004.10.008.
  * Garcia et al. [2019] D. Garcia, D. Pardo, V. M. Calo, Refined isogeometric analysis for fluid mechanics and electromagnetics, Computer Methods in Applied Mechanics and Engineering 356 (2019) 598–628, doi:10.1016/j.cma.2019.06.011.
  * Paszyński et al. [2018] M. Paszyński, L. Siwik, M. Woźniak, Concurrency of three-dimensional refined isogeometric analysis, Parallel Computing 80 (2018) 1–22, doi:10.1016/j.parco.2018.10.005.
  * Siwik et al. [2019] L. Siwik, M. Woźniak, V. Trujillo, D. Pardo, V. M. Calo, M. Paszyński, Parallel Refined Isogeometric Analysis in 3D, IEEE Transactions on Parallel and Distributed Systems 30 (5) (2019) 1134–1142, doi:10.1109/tpds.2018.2879664.
  * Hashemian et al. [2021a] A. Hashemian, D. Pardo, V. M. Calo, Refined isogeometric analysis for generalized Hermitian eigenproblems, Computer Methods in Applied Mechanics and Engineering 381 (2021a) 113823, doi:10.1016/j.cma.2021.113823.
  * Hashemian et al. [2021b] A. Hashemian, D. Garcia, J. A. Rivera, D. Pardo, Massive database generation for 2.5D borehole electromagnetic measurements using refined isogeometric analysis, Computers & Geosciences 155 (2021b) 104808, doi:10.1016/j.cageo.2021.104808.
  * Cottrell et al. [2006] J. A. Cottrell, A. Reali, Y. Bazilevs, T. J. R. Hughes, Isogeometric analysis of structural vibrations, Computer Methods in Applied Mechanics and Engineering 195 (49-50) (2006) 5257–5296, doi:10.1016/j.cma.2005.09.027.
  * Hughes et al. [2008] T. J. R. Hughes, A. Reali, G. Sangalli, Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: Comparison of $p$-method finite elements with $k$-method NURBS, Computer Methods in Applied Mechanics and Engineering 197 (2008) 4104–4124, doi:10.1016/j.cma.2008.04.006.
  * Buffa et al. [2010a] A. Buffa, G. Sangalli, R. Vázquez, Isogeometric analysis in electromagnetics: B-splines approximation, Computer Methods in Applied Mechanics and Engineering 199 (17-20) (2010a) 1143–1152, doi:10.1016/j.cma.2009.12.002.
  * Hughes et al. [2014] T. J. R. Hughes, J. A. Evans, A. Reali, Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems, Computer Methods in Applied Mechanics and Engineering 272 (2014) 290–320, doi:10.1016/j.cma.2013.11.012.
  * Puzyrev et al. [2017] V. Puzyrev, Q. Deng, V. Calo, Dispersion-optimized quadrature rules for isogeometric analysis: Modified inner products, their dispersion properties, and optimally blended schemes, Computer Methods in Applied Mechanics and Engineering 320 (2017) 421–443, doi:10.1016/j.cma.2017.03.029.
  * Hosseini et al. [2018] S. F. Hosseini, A. Hashemian, A. Reali, On the application of curve reparameterization in isogeometric vibration analysis of free-from curved beams, Computers & Structures 209 (2018) 117–129, doi:10.1016/j.compstruc.2018.08.009.
  * Mazza et al. [2019] M. Mazza, C. Manni, A. Ratnani, S. Serra-Capizzano, H. Speleers, Isogeometric analysis for 2D and 3D curl–div problems: Spectral symbols and fast iterative solvers, Computer Methods in Applied Mechanics and Engineering 344 (2019) 970–997, doi:10.1016/j.cma.2018.10.008.
  * Deng et al. [2019a] Q. Deng, V. Puzyrev, V. Calo, Optimal spectral approximation of $2n$-order differential operators by mixed isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 343 (2019a) 297–313, doi:10.1016/j.cma.2018.08.042.
  * Deng et al. [2019b] Q. Deng, V. Puzyrev, V. Calo, Isogeometric spectral approximation for elliptic differential operators, Journal of Computational Science 36 (2019b) 100879, doi:10.1016/j.jocs.2018.05.009.
  * Gao et al. [2020] H. Gao, L. Chen, H. Lian, C. Zheng, H. Xu, T. Matsumoto, Band structure analysis for 2D acoustic phononic structure using isogeometric boundary element method, Advances in Engineering Software 149 (2020) 102888, doi:10.1016/j.advengsoft.2020.102888.
  * Puzyrev et al. [2018] V. Puzyrev, Q. Deng, V. Calo, Spectral approximation properties of isogeometric analysis with variable continuity, Computer Methods in Applied Mechanics and Engineering 334 (2018) 22–39, doi:10.1016/j.cma.2018.01.042.
  * Sadler and Yang [1997] J. P. Sadler, Z. Yang, A comprehensive study of modal characteristics of a cylindrical manipulator with both link and joint flexibility, Mechanism and Machine Theory 32 (8) (1997) 941–956, doi:10.1016/S0094-114X(97)00001-3.
  * Olson and Vandini [1989] L. Olson, T. Vandini, Eigenproblems from finite element analysis of fluid–structure interactions, Computers & Structures 33 (3) (1989) 679–687, doi:10.1016/0045-7949(89)90242-3.
  * Sima et al. [2004] D. M. Sima, S. V. Huffel, G. H. Golub, Regularized Total Least Squares Based on Quadratic Eigenvalue Problem Solvers, BIT Numerical Mathematics 44 (4) (2004) 793–812, doi:10.1007/s10543-004-6024-8.
  * Cooke and Levush [2000] S. J. Cooke, B. Levush, Eigenmode Solution of 2-D and 3-D Electromagnetic Cavities Containing Absorbing Materials Using the Jacobi–Davidson Algorithm, Journal of Computational Physics 157 (1) (2000) 350–370, doi:10.1006/jcph.1999.6386.
  * Bermúdez et al. [2000] A. Bermúdez, R. G. Durán, R. Rodríguez, J. Solomin, Finite Element Analysis of a Quadratic Eigenvalue Problem Arising in Dissipative Acoustics, SIAM Journal on Numerical Analysis 38 (1) (2000) 267–291, doi:10.1137/s0036142999360160.
  * Campos and Roman [2016] C. Campos, J. E. Roman, Restarted Q-Arnoldi-type methods exploiting symmetry in quadratic eigenvalue problems, BIT Numerical Mathematics 56 (4) (2016) 1213–1236, doi:10.1007/s10543-016-0601-5.
  * Campos and Roman [2020] C. Campos, J. E. Roman, Inertia-based spectrum slicing for symmetric quadratic eigenvalue problems, Numerical Linear Algebra with Applications 27 (4) (2020) e2293, doi:10.1002/nla.2293.
  * Bai et al. [2000] Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, H. van der Vorst, Templates for the Solution of Algebraic Eigenvalue Problems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000\.
  * Tisseur and Meerbergen [2001] F. Tisseur, K. Meerbergen, The Quadratic Eigenvalue Problem, SIAM Review 43 (2) (2001) 235–286, doi:10.1137/s0036144500381988.
  * Vázquez and Buffa [2010] R. Vázquez, A. Buffa, Isogeometric Analysis for Electromagnetic Problems, IEEE Transactions on Magnetics 46 (8) (2010) 3305–3308, doi:10.1109/tmag.2010.2044563.
  * Evans and Hughes [2013] J. A. Evans, T. J. R. Hughes, Isogeometric divergence-conforming B-splines for the steady Navier–Stokes equations, Mathematical Models and Methods in Applied Sciences 23 (08) (2013) 1421–1478, doi:10.1142/s0218202513500139.
  * Buffa et al. [2010b] A. Buffa, C. de Falco, G. Sangalli, IsoGeometric Analysis: Stable elements for the 2D Stokes equation, International Journal for Numerical Methods in Fluids 65 (11-12) (2010b) 1407–1422, doi:10.1002/fld.2337.
  * Qian and Lin [2007] J. Qian, W.-W. Lin, A numerical method for quadratic eigenvalue problems of gyroscopic systems, Journal of Sound and Vibration 306 (1-2) (2007) 284–296, doi:10.1016/j.jsv.2007.05.009.
  * Zhang [2018] Z. Zhang, Comparison results for eigenvalues of curl curl operator and Stokes operator, Zeitschrift für angewandte Mathematik und Physik 69 (4) (2018) 104, doi:10.1007/s00033-018-0997-7.
  * Piegl and Tiller [1997] L. Piegl, W. Tiller, The NURBS Book, Springer-Verlag, New York, NY, 2nd edn., 1997.
  * Demkowicz et al. [2000] L. Demkowicz, P. Monk, L. Vardapetyan, W. Rachowicz, De Rham diagram for $hp$ finite element spaces, Computers & Mathematics with Applications 39 (7-8) (2000) 29–38, doi:10.1016/s0898-1221(00)00062-6.
  * Karypis and Kumar [1998] G. Karypis, V. Kumar, A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs, SIAM Journal on Scientific Computing 20 (1) (1998) 359–392, doi:10.1137/s1064827595287997.
  * Ericsson and Ruhe [1980] T. Ericsson, A. Ruhe, The Spectral Transformation Lanczos Method for the Numerical Solution of Large Sparse Generalized Symmetric Eigenvalue Problems, Mathematics of Computation 35 (152) (1980) 1251–1268, doi:10.2307/2006390.
  * Nour-Omid et al. [1987] B. Nour-Omid, B. N. Parlett, T. Ericsson, P. S. Jensen, How to implement the spectral transformation, Mathematics of Computation 48 (178) (1987) 663–673, doi:10.1090/S0025-5718-1987-0878698-5.
  * Grimes et al. [1994] R. G. Grimes, J. G. Lewis, H. D. Simon, A Shifted Block Lanczos Algorithm for Solving Sparse Symmetric Generalized Eigenproblems, SIAM Journal on Matrix Analysis and Applications 15 (1) (1994) 228–272, doi:10.1137/S0895479888151111.
  * Xue and Elman [2011] F. Xue, H. C. Elman, Fast inexact subspace iteration for generalized eigenvalue problems with spectral transformation, Linear Algebra and its Applications 435 (3) (2011) 601–622, doi:10.1016/j.laa.2010.06.021.
  * Stewart [2001] G. W. Stewart, Matrix Algorithms, Volume II: Eigensystems, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2001.
  * Wu and Simon [2000] K. Wu, H. Simon, Thick-Restart Lanczos Method for Large Symmetric Eigenvalue Problems, SIAM Journal on Matrix Analysis and Applications 22 (2) (2000) 602–616, doi:10.1137/S0895479898334605.
  * Stewart [2002a] G. W. Stewart, A Krylov–Schur Algorithm for Large Eigenproblems, SIAM Journal on Matrix Analysis and Applications 23 (3) (2002a) 601–614, doi:10.1137/S0895479800371529.
  * Stewart [2002b] G. W. Stewart, Addendum to “A Krylov–Schur Algorithm for Large Eigenproblems”, SIAM Journal on Matrix Analysis and Applications 24 (2) (2002b) 599–601, doi:10.1137/S0895479802403150.
  * Campos and Roman [2012] C. Campos, J. E. Roman, Strategies for spectrum slicing based on restarted Lanczos methods, Numerical Algorithms 60 (2012) 279–295, doi:10.1007/s11075-012-9564-z.
  * Parlett [1998] B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice-Hall, Inc., Upper Saddle River, NJ, 1998.
  * Nakatsukasa and Noferini [2019] Y. Nakatsukasa, V. Noferini, Inertia laws and localization of real eigenvalues for generalized indefinite eigenvalue problems, Linear Algebra and its Applications 578 (2019) 272–296, doi:10.1016/j.laa.2019.05.010.
  * Guo and Lancaster [2005] C.-H. Guo, P. Lancaster, Algorithms for hyperbolic quadratic eigenvalue problems, Mathematics of Computation 74 (252) (2005) 1777–1792, doi:10.1090/s0025-5718-05-01748-5.
  * Niendorf and Voss [2010] V. Niendorf, H. Voss, Detecting hyperbolic and definite matrix polynomials, Linear Algebra and its Applications 432 (4) (2010) 1017–1035, doi:10.1016/j.laa.2009.10.014.
  * Sleijpen et al. [1996] G. L. G. Sleijpen, A. G. L. Booten, D. R. Fokkema, H. A. van der Vorst, Jacobi-davidson type methods for generalized eigenproblems and polynomial eigenproblems, BIT Numerical Mathematics 36 (3) (1996) 595–633, doi:10.1007/bf01731936.
  * Bai and Su [2005] Z. Bai, Y. Su, SOAR: A Second-order Arnoldi Method for the Solution of the Quadratic Eigenvalue Problem, SIAM Journal on Matrix Analysis and Applications 26 (3) (2005) 640–659, doi:10.1137/s0895479803438523.
  * Lu et al. [2016] D. Lu, Y. Su, Z. Bai, Stability Analysis of the Two-level Orthogonal Arnoldi Procedure, SIAM Journal on Matrix Analysis and Applications 37 (1) (2016) 195–214, doi:10.1137/151005142.
  * Collier et al. [2013] N. Collier, L. Dalcin, D. Pardo, V. M. Calo, The Cost of Continuity: Performance of Iterative Solvers on Isogeometric Finite Elements, SIAM Journal on Scientific Computing 35 (2) (2013) A767–A784, doi:10.1137/120881038.
  * Sarmiento et al. [2017] A. F. Sarmiento, A. M. A. Côrtes, D. A. Garcia, L. Dalcin, N. Collier, V. M. Calo, PetIGA-MF: A multi-field high-performance toolbox for structure-preserving B-splines spaces, Journal of Computational Science 18 (2017) 117–131, doi:10.1016/j.jocs.2016.09.010.
  * Dalcin et al. [2016] L. Dalcin, N. Collier, P. Vignal, A. M. A. Côrtes, V. M. Calo, PetIGA: A framework for high-performance isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 308 (2016) 151–181, doi:10.1016/j.cma.2016.05.011.
  * Balay et al. [1997] S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith, Efficient Management of Parallelism in Object Oriented Numerical Software Libraries, in: E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkhäuser Press, 163–202, 1997\.
  * Vignal et al. [2015] P. Vignal, A. Sarmiento, A. M. A. Côrtes, L. Dalcin, V. M. Calo, Coupling Navier-Stokes and Cahn-Hilliard Equations in a Two-dimensional Annular flow Configuration, Procedia Computer Science 51 (2015) 934–943, doi:10.1016/j.procs.2015.05.228.
  * Espath et al. [2016] L. F. R. Espath, A. F. Sarmiento, P. Vignal, B. O. N. Varga, A. M. A. Cortes, L. Dalcin, V. M. Calo, Energy exchange analysis in droplet dynamics via the Navier–Stokes–Cahn–Hilliard model, Journal of Fluid Mechanics 797 (2016) 389–430, doi:10.1017/jfm.2016.277.
  * Côrtes et al. [2017] A. M. A. Côrtes, L. Dalcin, A. F. Sarmiento, N. Collier, V. M. Calo, A scalable block-preconditioning strategy for divergence-conforming B-spline discretizations of the Stokes problem, Computer Methods in Applied Mechanics and Engineering 316 (2017) 839–858, doi:10.1016/j.cma.2016.10.014.
  * Espath et al. [2017] L. F. R. Espath, A. F. Sarmiento, L. Dalcin, V. M. Calo, On the thermodynamics of the Swift–Hohenberg theory, Continuum Mechanics and Thermodynamics 29 (6) (2017) 1335–1345, doi:10.1007/s00161-017-0581-y.
  * Sarmiento et al. [2018] A. F. Sarmiento, L. F. R. Espath, P. Vignal, L. Dalcin, M. Parsani, V. M. Calo, An energy-stable generalized-$\alpha$ method for the Swift–Hohenberg equation, Journal of Computational and Applied Mathematics 344 (2018) 836–851, doi:10.1016/j.cam.2017.11.004.
  * Clavijo et al. [2019] S. P. Clavijo, A. F. Sarmiento, L. F. R. Espath, L. Dalcin, A. M. A. Cortes, V. M. Calo, Reactive $n$-species Cahn–Hilliard system: A thermodynamically-consistent model for reversible chemical reactions, Journal of Computational and Applied Mathematics 350 (2019) 143–154, doi:10.1016/j.cam.2018.10.007.
  * Hernandez et al. [2005] V. Hernandez, J. E. Roman, V. Vidal, SLEPc: A Scalable and Flexible Toolkit for the Solution of Eigenvalue Problems, ACM Transactions on Mathematical Software 31 (3) (2005) 351–362, doi:10.1145/1089014.1089019.
  * Romero and Roman [2014] E. Romero, J. E. Roman, A parallel implementation of Davidson methods for large-scale eigenvalue problems in SLEPc, ACM Transactions on Mathematical Software 40 (2) (2014) 1–29, doi:10.1145/2543696.
  * Faber et al. [2018] B. J. Faber, M. J. Pueschel, P. W. Terry, C. C. Hegna, J. E. Roman, Stellarator microinstabilities and turbulence at low magnetic shear, Journal of Plasma Physics 84 (5) (2018) 905840503, doi:10.1017/s0022377818001022.
  * Keçeli et al. [2018] M. Keçeli, F. Corsetti, C. Campos, J. E. Roman, H. Zhang, Á. Vázquez-Mayagoitia, P. Zapol, A. F. Wagner, SIESTA-SIPs: Massively parallel spectrum-slicing eigensolver for an ab initio molecular dynamics package, Journal of Computational Chemistry 39 (22) (2018) 1806–1814, doi:10.1002/jcc.25350.
  * Araujo C. et al. [2020] J. C. Araujo C., C. Campos, C. Engström, J. E. Roman, Computation of scattering resonances in absorptive and dispersive media with applications to metal-dielectric nano-structures, Journal of Computational Physics 407 (2020) 109220, doi:10.1016/j.jcp.2019.109220.
  * Amestoy et al. [2001] P. R. Amestoy, I. S. Duff, J.-Y. L’Excellent, J. Koster, A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling, SIAM Journal on Matrix Analysis and Applications 23 (1) (2001) 15–41, doi:10.1137/s0895479899358194.
  * Shahriari et al. [2021] M. Shahriari, D. Pardo, J. A. Rivera, C. Torres-Verdín, A. Picon, J. Del Ser, S. Ossandón, V. M. Calo, Error control and loss functions for the deep learning inversion of borehole resistivity measurements, International Journal for Numerical Methods in Engineering 122 (6) (2021) 1629–1657, doi:10.1002/nme.6593.
  * Nicolet and Geuzaine [2006] A. Nicolet, C. Geuzaine, Waveguide Propagation Modes and Quadratic Eigenvalue Problems, in: Proceedings of 6th International Conference on Computational Electromagnetics, IEEE, Frankfurt am Main, Germany, 2006\.
  * Bermúdez and Rodríguez [1999] A. Bermúdez, R. Rodríguez, Modelling and numerical solution of elastoacoustic vibrations with interface damping, International Journal for Numerical Methods in Engineering 46 (10) (1999) 1763–1779, doi:10.1002/(sici)1097-0207(19991210)46:10<1763::aid-nme723>3.0.co;2-6.
  * Tisseur [2000] F. Tisseur, Backward error and condition of polynomial eigenvalue problems, Linear Algebra and its Applications 309 (1) (2000) 339–361, doi:10.1016/S0024-3795(99)00063-4.