Research

astrophysics, cosmology, time-domain astronomy

Ariel Goobar

# Strongly lensed supernovae: lessons learned

A.Goobar  J. Johansson and A. Sagués Carracedo  The Oskar Klein Centre,
Department of Physics, Stockholm University, SE 106 91 Stockholm, Sweden
<EMAIL_ADDRESS>
###### Abstract

Since a few years, we have finally entered the era of discoveries of multiply-
imaged gravitationally lensed supernovae. To date, all cluster lensed
supernovae have been found from space, while those deflected by individual
galaxies were identified with wide-field ground-based surveys through the
magnification of "standard candles" method, i.e., without the need of
spatially resolving the individual images. We review the challenges in
identifying these extremely rare events, as well as the unique opportunities
they offer for time-delay cosmography and the study of the properties of the
deflecting bodies acting as lenses.

###### keywords:

supernovae, gravitational lensing, time-domain surveys

##  1 Introduction

Following the early days of supernova cosmology, where CCD cameras exceeding a
few arcminutes scale became available to search for high-redshift supernovae,
targeted observations of massive lensing clusters used as “gravitational
telescopes” were proposed to search for the most distant supernovae [1, 2, 3].
The gain factor in exposure length is $\mu^{2}$, where $\mu$ is the flux
amplification provided by the lens. However, this is partially balanced as the
solid angle at the source planes shrinks by a factor $\mu$ behind the lens.
Hence, earlier attempts using ground-based optical and infrared instruments to
do monthly cadenced observations of lensing clusters did not uncover any
multiply-imaged supernova [4, 5]. It was only through Hubble Space Telescope
(HST) observations that the first cluster lensed supernova, a core-collapse
supernova at $z=1.49$ lensed by the MACS J1149.6+2223 cluster was found [6].
The supernova was named "SN Refsdal", honouring the memory of Sjur Refsdal who
first proposed to use time delays between multiple images of strongly lensed
supernova (gLSNe) to measure the Hubble constant [7]. Continued HST monitoring
of massive clusters, and more recently also with JWST, has led to the
discovery of several other multiply imaged supernovae behind clusters [8, 9,
10, 11], including three Type Ia supernovae (SNe Ia): SN H0pe and two
"siblings", i.e., SNe hosted by the same galaxy, SN Requiem and SN Encore.
While cluster lensed SNe have so far only been discovered from space, the
three galaxy lensed SNe found to date were found in very wide-field surveys
with ground-based telescopes. PS1-10afx was first reported as an unusual
superluminous SN in the PanSTARRS transient survey [12]. Three years after the
SN discovery (and too late for high-spatial resolution follow-up), it was
shown in [13] that it was a highly magnified SN Ia at redshift $z=1.388$, and
eventually also the lens at $z=1.117$ was identified [14]. Since then, two
multiply-imaged SNe Ia have been found at Palomar Observatory. Starting with
iPTF16geu, a SNIa at redshift $z_{s}=0.409$, deflected by a galaxy at
$z_{l}=0.2163$, detected by the intermediate Palomar Transient Factory [15]
and followed by another SNIa, SN Zwicky ($z_{s}=0.354;z_{l}=0.226$) [16] by
the ongoing Zwicky Transient Facility (ZTF). Figure 1 shows space imaging for
three multiply-imaged SNe Ia, highlighting the different angular scales for
cluster and galaxy lens systems found to date. For a recent review of the
status of lensing of supernovae, see [17]. The current manuscript focuses on
lessons learned on gLSNe findings from ground-based transient surveys,
iPTF16geu and SN Zwicky in particular.

Figure 1: Gravitationally lensed Type Ia supernovae with multiple images, SN
H0pe [10] lensed by a cluster of galaxies and iPTF16geu [15] and SN Zwicky
[16, 18] by individual galaxies. For the latter, the image flux ratios suggest
that additional micro- or millilensing from stellar objects or substructures
is taking place in the deflecting galaxy [19, 20]. For iPTF16geu, a
significant part of the intensity differences are due to extinction in the
lensing galaxy [21].

##  2 Spatially unresolved strongly lensed SNe

Wide-field imaging transient surveys like Palomar’s PTF (2009-2012), iPTF
(2013-2017) and ZTF (operating since 2018) have the ability to cover the
entire visible sky from the Northern hemisphere in a single night. The
extremely large search area of ZTF, facilitated by its 47.sq. deg field-of-
view camera, makes it especially suitable to detect rare transient phenomena.
The limiting factors are the small collecting area of the 1.2m telescope and
the very coarse spatial resolution. With 1$"$ pixel plate-scale and typically
2$"$ seeing at Palomar, detecting spatially resolved gLSNe would be extremely
rare. Simulations of the ZTF survey [22] indicate that only about 2% of the
gLSNe within discovery range from ZTF would have image separations exceeding
3$"$, at which point they could be detected as two individual point sources.
The time delay between images is typically shorter than the typical time scale
of the lightcurves, hence making also quite challenging to identify gLSNe from
the vast pool of regular SN lightcurves through multiple detections separated
in time. Figure 2 shows the distribution of time delays between SN images and
the characteristic angular scale of strong lensing, the Einstein radius
$\theta_{E}$ expected from simulations of the ZTF survey [23]. The extremely
compact multi-image systems iPTF16geu $\theta_{E}=0.3"$ and SN Zwicky
$\theta_{E}=0.16"$ (shown in Figure 1) could be identified through the
magnification method. As indicated with stars in Figure 2, these two SNe where
highly magnified $\Delta m=2.5\log_{10}(\mu)>3$mag, split into four images.
Thanks to the "standard candle" nature of Type Ia supernovae, typically
showing some $\sim 0.15$ mag of scatter after applying corrections for
lightcurve stretch and colour, can be easily identified as outliers in
brightness, provided a spectrum is available with the spectral classification
and the redshift of the SN. We will return to this issue in Section 4.

Figure 2: Left: Probability distribution of the time delays (in days) between
multiple SN images vs the Einstein radius (arcseconds) for gLSN systems
expected in ZTF. Right: Magnification (in magnitudes) vs Einstein radius. The
two stars show the location in the parameter space of iPTF16geu and SN Zwicky.
The black contours indicate the 68% and 95% confidence regions.

##  3 A different population of lens systems

As can be appreciated from Figures 1 and 2, the systems found to date form the
ground are extremely compact. Such small angular separation lensing systems
are rarely found by other means given the extreme spatial resolution needed.
Hence, it was shown in [16] that the gLSNe uncover a yet unexplored population
of low stellar mass lensing galaxies. In particular, the compact systems
provide interesting insights into the inner $\sim$1 kpc region of lensing
galaxies. The downside, discussed further in the Section 5, is that they are
not suitable for cosmographic time-delay measurements.

##  4 The discovery bottleneck: spectroscopic follow-up

Because of limited spectroscopic resources, only a small fraction of the
transient discoveries in iPTF and ZTF were followed-up with the necessary
spectroscopic screening needed to identify a lensed SN. While the photometric
detection threshold in ZTF is around 20.5-21 mag [24], the spectroscopic
classification, as a part of the Bright Transient Survey (BTS) is only
complete to 18.5 mag [25]. Simulations of the ZTF survey [23] (see also [26])
show that the bright threshold of the BTS spectroscopic classification has
been the bottleneck for identifying the gLSNe with the magnification method.
In a recent work [27], an archival study of ZTF data was used to search for
possible missed "live" candidates due to the magnitude limitation from BTS.
The search efficiency was enhanced by having access to galaxy redshifts from
the Dark Energy Spectroscopic Instrument (DESI), spatially associated with
longlived, red candidates. The search has produced a handful of intriguing
candidates. While superluminous supernovae (SLSNe) cannot be fully rejected as
a possible explanation, two archival ZTF events, are significantly different
from typical SLSNe and their lightcurves can be modelled as two-image lensed
SNIa systems. From this two-image modelling, time delays of $22\pm 3$ and
$34\pm 1$ days were estimated, respectively. If confirmed, it suggests that we
may have found the first events with longer time delays with ground-based
resources! The findings are in good agreement with the rate expectations from
survey simulations in [23], shown in Figure 3.

Figure 3: Expected yearly discovery rate as a function of the apparent
magnitude threshold for SNe Ia (blue) and CCSNe (red) lensed supernovae. The
dashed fainter curves are for subtypes of that are individual components of
the red curve for CC: IIP, IIn, and Ibc, of which IIn is the dominant type
wíth standard assumptions on their luminosity function and constant fraction
of CC population, independent of redshift. The vertical dashed lines indicate
the magnitude cuts of 18.5, 19, and 20 mag with the corresponding percentage
of each SN type up to that cut. 18.5 mag corresponds to the BTS magnitude
completeness cut, under which most supernovae are spectroscopically
classified. We also include a region up to 19 mag as BTS extends to such
magnitudes when the schedule allows is. Adapted from [23].

##  5 Time delays and the quest for the Hubble constant

One of the main motivations behind searching for gLSNe is to use their
lightcurves to measure the time delays between the multiple images, from which
the Hubble constant (H0) can be inferred, as first suggested by Refsdal in
1964 [7]. In recent years, the interest in this type of measurement has gained
a lot of interest due to the emergence of the so called "Hubble tension",
suggesting that the value of H0 obtained from the early universe CMB
anisotropy data, extrapolated to the present universe using the $\Lambda$CDM
model ($67.4\pm 0.5$ km s-1Mpc-1), is in conflict with the local distance
ladder measurement from the SH0ES team ($73.0\pm 1.0$ km s-1Mpc-1), see [28]
for a recent review. Time-delay cosmography offers an interesting independent
way to measure the Hubble constant and could provide further support or reject
the notion that physics beyond the $\Lambda$CDM model is required. For many
years, time-delay cosmography has been carried out exclusively with multiply-
imaged quasars, but the results are as of yet inconclusive (see e.g. [29] for
a status update). The smooth lightcurves of supernovae coupled with their
favourable time scales make them potentially superior to QSOs for time-delay
cosmography. Furthermore, unlike QSOs, supernovae fade in roughly a year time-
scale, allowing for detailed studies of the lens without contamination from
the lensed images. Hence, the possibility to complement the time-delay
cosmography from QSOs with gLSNe has generated a lot of interest.

Furthermore, thanks to the standard-candle nature of SNe Ia (after corrections
for colour and lightcurve shape), their magnification can be inferred up to an
uncertainty related to their intrinsic luminosity scatter, about 0.15 mag.
This is potentially a key feature, since it can be used to break the so called
mass-sheet degeneracy. In brief, the presence of a constant sheet of surface
mass density leaves the predicted images unchanged, but alters the time delay
between the images [30]. Breaking the mass-sheet degeneracy, e.g., through the
model independent measurement of the SNIa magnification, is therefore a very
important element for constraining H0 [31].

In the following sections we will discuss some additional challenges to break
the mass-sheet degeneracy posed by extinction by dust in the host and
deflecting galaxy, as well as micro and millilensing.

###  5.1 Cluster scale lenses

Through the monitoring of the multiple images of SN Refsdal [32] the time
delays and magnification ratios among the images were measured, and the most
accurate time delay between a pair of images was $376.0^{+5.6}_{-5.5}$ days
[33]. This time-delay measurement with a relative uncertainty of 1.5% provided
the first precise $H_{0}$ measurement from lensed SNe. Through lensing models
of the cluster, [33] found $H_{0}=66.6^{+4.1}_{-3.3}$ km s-1Mpc-1. More
recently, H0 was measured for SN H0pe. A combination of a spectroscopic [34]
(see Section 5.4) and photometric [35] time-delay measurement were compared to
the predictions of many cluster lens models to measure a value for the Hubble
constant [36]. In combination with the magnification of this SNIa, yielded a
value of $H_{0}=75.4^{+8.1}_{-5.5}$ km s-1Mpc-1.

###  5.2 Galaxy scale lenses

While the H0 measurements from SN Refsdal and SN H0pe are very encouraging and
exciting, some caution needs to be exercised in interpreting these results.
The lensing models are very challenging since the mass distributions of
clusters are quite complex, implying that the multiple image region of
clusters is expected to be rich in substructures. For that reason, galaxy
lenses are much simpler to model and therefore preferable, as they involve
smaller systematic uncertainties. However, the highly magnified compact
systems within reach for shallow surveys like ZTF are expected to produce
images which can be separated by just a few days, as shown in Figure 4, making
it rather challenging for precision measurements of time delays. That was the
case for both iPTF16geu [21] and SN Zwicky [16], as outlined in the next
Section.

Figure 4: Expected cumulative distribution of time delays for gLSNe discovered
by ZTF showing that the median time delay is close to 10 days for both core-
collapse and Type Ia supernovae. Adapted from [23].

###  5.3 The second maximum in the SN Ia near-IR lightcurves

Besides their "standard candle" nature, SNe Ia offer other benefits for time-
delay cosmography. For the restframe lightcurves in bands beyond the $r$
filter, a secondary maximum within about a month from the restframe $B$-band
lightcurve peak can be used to measure photometric time delays. This is
extremely useful, as it means that accurate arrival time differences between
SN images can be measured accurately, even when the first maximum is missed or
poorly sampled, as was the case for both iPTF16geu and SN H0pe, shown in
Figure 5.

Figure 5: Top row: Photometric time-delay estimates from resolved images of SN
H0pe (left) and iPTF16geu (right). In both cases, the "second IR maximum" was
used to measure the difference in arrival times between the multiple SN
images. Bottom row: spectroscopic time delays where the spectral features are
dated using the SED template in [37]. Observations and further details can be
found in [35, 34, 21, 19].

###  5.4 Spectroscopic time delays

The two bottom panels of Figure 5 show another unique feature of the use of
supernovae for time-delay cosmography. During the early phases of the
supernova explosion, the atmosphere is thick and the supernova spectral energy
distribution is formed by the outer layers with lower opacity, the
photosphere. For a homologous expansion, $r=v\cdot t$, that corresponds to
very high velocities. As the expansion thins out the atmosphere, the
photosphere recedes, and the typical absorption features come from closer the
centre, hence lower velocities. This change of velocities of the SN features
can be used to extract the phase of the supernova at the time of observations.
Spectroscopic time-delay measurements have been carried out successfully for
iPTF16geu [38] and SN H0pe [34], as shown in Figure 5.

###  5.5 Measuring time delays with (mostly) unresolved data

One of the important recent developments is the realisation that time-delays
can be inferred from the unresolved lightcurves, provided there is at least
one high-spatial resolution image of the system that gives the image
multiplicity, their positions, and the relative image fluxes, as shown in
Figure 6 from [16]. The publicly available, python-based software sntd [39]
was used for for inferring the restframe $B$-peak magnitude, the lightcurve
shape and colour SNIa SALT2 parameters [40] and the time-delays between the
images. Unresolved photometry from the Palomar and Liverpool telescopes in
$g,r,i,z$ filters were included in the fit, along with a model including the
flux contributions from the four sets of lightcurves accounting for
extinction, each one with their own time of maximum The fit is constrained by
imposing a prior on the image ratios at the date of the Keck/NIRC2
observations, shown in the right-hand side panel of Figure 6. The total
lensing magnification was fitted, $\mu=24.3\pm 2.7$. Although negligible time
delays were found in this case, the method is very promising for future
systems.

Figure 6: Ground-based (unresolved) lightcurves for SN Zwicky and resolved
Keck/AO images from [16]. Using the prior from the Keck/AO image showing the
quad configuartion and the fluxes at one epoch, time delays could be measured
accurately with the unresolved multi-band data.

##  6 The interstellar medium and differential extinction

As the light of gLSNe pierces through the inner regions of the deflecting
galaxies, measurements of magnification crucially depend on the ability to
accurately correct for losses due to scattering on dust grains, both in the
host and lensing galaxy. For that purpose, multi-band imaging is used, since
the magnitude increase (i.e., loss in flux in logarithmic units) to dust
extinction is roughly inversely proportional with wavelength,
$A_{\lambda}\propto\lambda^{-1}$. To further complicate matters, observational
evidence suggests that the composition and grain size distribution in the
extragalactic interstellar medium could be very diverse [41], with total-to-
selective extinction $R_{V}=A_{V}/E(B-V)$ potentially quite different from the
Milky-Way value, hence the need to both fit the colour excess $E(B-V)$ and
$R_{V}$, even from individual images. Since precise information of the range
of properties for dimming by dust in other galaxies is so critical for
accurate distance measurements in cosmology, it is very exciting to be able to
carry out such measurements with resolved images of lensed SNe, as was the
case for iPTF16geu [21], shown in Figure 7. Images C and D (see Figure 1)
where particularly reddened, and thanks to having HST images in at least four
filters useful constraints could be set on both the colour excess and $R_{V}$,
providing unique a test of dust grain density and properties multiple lines of
sight in an intermediate redshift galaxy.

Figure 7: Inferred wavelength dependent extinction, $A_{\lambda}$, for the
four resolved images of iPTF16geu measured with HST. The absorption from the
host galaxy dust is plotted with dotted black line. For Image A we can see
that the host galaxy is the dominant source of extinction, and for Images B,
C, D there is a progressively larger contribution from the dust in the lens
galaxy (see [21] for data and analysis information).

##  7 Macro vs milli/microlensing

Using multi-band follow-up observations of iPTF16geu with HST, an accurate
(model independent) measurement of the total magnification was made,
$\mu=67.8^{+2.6}_{-2.9}$ [21], after correction for non-negligible extinction
by dust in both the host and lens galaxies, as discussed in Section 6. The
time delays between the SN images for this system were very small, about a day
or less [42, 21]. However, the flux ratios between the supernova images (see
Fig. 1) were not consistent with expectations from (macro) lensing of a smooth
extended deflector, even accounting for differential extinction, hinting at
additional lensing contributions from galactic sub halos (millilensing) or
stellar objects (microlensing) [19]. In either case, the sub-spliting of the
SN images is too small to be resolved as it is of order milli-arcseconds or
less, whereas e.g., the HST angular resolution is at least a few tens of
milli-arcseconds. The situation was very similar for SN Zwicky ($\mu=23.7\pm
3.2$), except that there was no ambiguity between dimming by dust and
microlensing (de-)magnification [16, 18]. While extinction by dust can be
identified and corrected through its wavelength dependence, microlensing by
stellar objects is a more severe challenge [43]. The intrinsic size of a SN is
comparable to the Einstein radius of an individual star in the deflecting
galaxy. Hence, the observed magnification is sensitive to the unknown
positions of stars and substructures in the lensing galaxy. A thorough
discussion on microlensing of SNe can be found in [17] and potential means to
mitigate this issue has been discussed in e.g., [44, 45]. On the positive
side, the image flux ratios observed can be used to infer limits on the
possible dark matter contribution from e.g., primordial black holes over a
wide mass range [20].

##  8 gLSNe in the LSST era

Thanks to the wide-field coverage and faint photometric limit of the Legacy
Survey of Space and Time (LSST) survey at the Vera Rubin Observatory, the
discovery rate of strongly lensed supernovae is expected to increase
dramatically. Simulation studies have shown that hundreds of lensed SNe should
be found [22], many of which will be spatially resolved [46]. For SNe Ia, [47]
found that a ‘gold sample’ of $\sim$ 10 lensed SNIa per year can be expected,
with time delay above 10 days caught before peak, and sufficiently bright
(below 22.5 mag) for spectroscopic follow-up observations. In three years of
LSST operations, such a sample can yield a 1.5% measurement of the Hubble
constant.

##  9 Conclusion

The rapid developments in time-domain astronomy, including wide-field imaging
from the ground and very sensitive near-IR space instruments has led to an
exciting development in the discovery of strongly lensed supernovae: the era
of time-delay cosmography with supernovae has begun! Ground-based searches
have uncovered a population of compact lens systems, where micro and
millilensing effects are very important to characterise, also since they
provide tests for the nature of dark matter. Lessons learned from searches to
date give rise to great optimism, as instruments to be deployed in the
immediate future will greatly enhance the feasibility and science reach of
cosmology and astrophysics with lensed supernovae.

Many thanks to the Royal Society for the great hospitality during the meeting.
The authors acknowledge the work behind these results from whole the ZTF
lensed supernova working group as well as financial support from
Vetenskapsrådet, the Swedish Research Council through grants project Dnr
2020-03444, the G.R.E.A.T research environment, Dnr 2016-06012, and the
Swedish National Space Agency, Dnr 2023-00226.

## References

  * [1] Kovner I, Paczynski B. 1988 Supernovae in Luminous Arcs. ApJ 335, L9. (10.1086/185328)
  * [2] Sullivan M, Ellis R, Nugent P, Smail I, Madau P. 2000 A strategy for finding gravitationally lensed distant supernovae. MNRAS 319, 549–556. (10.1046/j.1365-8711.2000.03875.x)
  * [3] Gunnarsson C, Goobar A. 2003 Massive galaxy clusters as gravitational telescopes for distant supernovae. A&A 405, 859–866. (10.1051/0004-6361:20030648)
  * [4] Goobar A, Paech K, Stanishev V, Amanullah R, Dahlén T, Jönsson J, Kneib JP, Lidman C, Limousin M, Mörtsell E, Nobili S, Richard J, Riehm T, von Strauss M. 2009 Near-IR search for lensed supernovae behind galaxy clusters. II. First detection and future prospects. A&A 507, 71–83. (10.1051/0004-6361/200811254)
  * [5] Amanullah R, Goobar A, Clément B, Cuby JG, Dahle H, Dahlén T, Hjorth J, Fabbro S, Jönsson J, Kneib JP, Lidman C, Limousin M, Milvang-Jensen B, Mörtsell E, Nordin J, Paech K, Richard J, Riehm T, Stanishev V, Watson D. 2011 A Highly Magnified Supernova at z = 1.703 behind the Massive Galaxy Cluster A1689. ApJ 742, L7. (10.1088/2041-8205/742/1/L7)
  * [6] Kelly PL, Rodney SA, Treu T, Foley RJ, Brammer G, Schmidt KB, Zitrin A, Sonnenfeld A, Strolger LG, Graur O, Filippenko AV, Jha SW, Riess AG, Bradac M, Weiner BJ, Scolnic D, Malkan MA, von der Linden A, Trenti M, Hjorth J, Gavazzi R, Fontana A, Merten JC, McCully C, Jones T, Postman M, Dressler A, Patel B, Cenko SB, Graham ML, Tucker BE. 2015 Multiple images of a highly magnified supernova formed by an early-type cluster galaxy lens. Science 347, 1123–1126. (10.1126/science.aaa3350)
  * [7] Refsdal S. 1964 On the possibility of determining Hubble’s parameter and the masses of galaxies from the gravitational lens effect. MNRAS 128, 307. (10.1093/mnras/128.4.307)
  * [8] Rodney SA, Brammer GB, Pierel JDR, Richard J, Toft S, O’Connor KF, Akhshik M, Whitaker K. 2021 A Gravitationally Lensed Supernova with an Observable Two-Decade Time Delay. arXiv e-prints p. arXiv:2106.08935.
  * [9] Chen W, Kelly PL, Oguri M, Broadhurst TJ, Diego JM, Emami N, Filippenko AV, Treu TL, Zitrin A. 2022 Shock cooling of a red-supergiant supernova at redshift 3 in lensed images. Nature 611, 256–259. (10.1038/s41586-022-05252-5)
  * [10] Frye BL, Pascale M, Pierel J, Chen W, Foo N, Leimbach R, Garuda N, Cohen S, Kamieneski P, Windhorst R, Koekemoer AM, Kelly P, Summers J, Engesser M, Liu D, Furtak L, Polletta M, Harrington K, Willner S, Diego JM, Jansen R, Coe D, Conselice CJ, Dai L, Dole H, D’Silva JCJ, Driver S, Grogin N, Marshall MA, Meena A, Nonino M, Ortiz, Rafael I, Pirzkal N, Robotham A, Ryan RE, Strolger L, Tompkins S, Trussler J, Willmer C, Yan H, Yun MS, Zitrin A. 2023 The JWST Discovery of the Triply-imaged Type Ia “Supernova H0pe” and Observations of the Galaxy Cluster PLCK G165.7+67.0. arXiv e-prints p. arXiv:2309.07326. (10.48550/arXiv.2309.07326)
  * [11] Pierel JDR, Newman AB, Dhawan S, Gu M, Joshi BA, Li T, Schuldt S, Strolger LG, Suyu SH, Caminha GB, Cohen SH, Diego JM, Dsilva JCJ, Ertl S, Frye BL, Granata G, Grillo C, Koekemoer AM, Li J, Robotham A, Summers J, Treu T, Windhorst RA, Zitrin A, Agarwal S, Agrawal A, Arendse N, Belli S, Burns C, Cañameras R, Chakrabarti S, Chen W, Collett TE, Coulter DA, Ellis RS, Engesser M, Foo N, Fox OD, Gall C, Garuda N, Gezari S, Gomez S, Glazebrook K, Hjorth J, Huang X, Jha SW, Kamieneski PS, Kelly P, Larison C, Moustakas LA, Pascale M, Pérez-Fournon I, Petrushevska T, Poidevin F, Rest A, Shahbandeh M, Shajib AJ, Siebert M, Storfer C, Talbot M, Wang Q, Wevers T, Zenati Y. 2024 Lensed Type Ia Supernova “Encore” at z=2: The First Instance of Two Multiply-Imaged Supernovae in the Same Host Galaxy. arXiv e-prints p. arXiv:2404.02139. (10.48550/arXiv.2404.02139)
  * [12] Chornock R, Berger E, Rest A, Milisavljevic D, Lunnan R, Foley RJ, Soderberg AM, Smartt SJ, Burgasser AJ, Challis P, Chomiuk L, Czekala I, Drout M, Fong W, Huber ME, Kirshner RP, Leibler C, McLeod B, Marion GH, Narayan G, Riess AG, Roth KC, Sanders NE, Scolnic D, Smith K, Stubbs CW, Tonry JL, Valenti S, Burgett WS, Chambers KC, Hodapp KW, Kaiser N, Kudritzki RP, Magnier EA, Price PA. 2013 PS1-10afx at z = 1.388: Pan-STARRS1 Discovery of a New Type of Superluminous Supernova. ApJ 767, 162. (10.1088/0004-637X/767/2/162)
  * [13] Quimby RM, Werner MC, Oguri M, More S, More A, Tanaka M, Nomoto K, Moriya TJ, Folatelli G, Maeda K, Bersten M. 2013 Extraordinary Magnification of the Ordinary Type Ia Supernova PS1-10afx. ApJ 768, L20. (10.1088/2041-8205/768/1/L20)
  * [14] Quimby RM, Oguri M, More A, More S, Moriya TJ, Werner MC, Tanaka M, Folatelli G, Bersten MC, Maeda K, Nomoto K. 2014 Detection of the Gravitational Lens Magnifying a Type Ia Supernova. Science 344, 396–399. (10.1126/science.1250903)
  * [15] Goobar A, Amanullah R, Kulkarni SR, Nugent PE, Johansson J, Steidel C, Law D, Mörtsell E, Quimby R, Blagorodnova N, Brandeker A, Cao Y, Cooray A, Ferretti R, Fremling C, Hangard L, Kasliwal M, Kupfer T, Lunnan R, Masci F, Miller AA, Nayyeri H, Neill JD, Ofek EO, Papadogiannakis S, Petrushevska T, Ravi V, Sollerman J, Sullivan M, Taddia F, Walters R, Wilson D, Yan L, Yaron O. 2017 iPTF16geu: A multiply imaged, gravitationally lensed type Ia supernova. Science 356, 291–295. (10.1126/science.aal2729)
  * [16] Goobar A, Johansson J, Schulze S, Arendse N, Carracedo AS, Dhawan S, Mörtsell E, Fremling C, Yan L, Perley D, Sollerman J, Joseph R, Hinds KR, Meynardie W, Andreoni I, Bellm E, Bloom J, Collett TE, Drake A, Graham M, Kasliwal M, Kulkarni SR, Lemon C, Miller AA, Neill JD, Nordin J, Pierel J, Richard J, Riddle R, Rigault M, Rusholme B, Sharma Y, Stein R, Stewart G, Townsend A, Vinko J, Wheeler JC, Wold A. 2023 Uncovering a population of gravitational lens galaxies with magnified standard candle SN Zwicky. Nature Astronomy 7, 1098–1107. (10.1038/s41550-023-01981-3)
  * [17] Suyu SH, Goobar A, Collett T, More A, Vernardos G. 2024 Strong Gravitational Lensing and Microlensing of Supernovae. Space Sci. Rev. 220, 13. (10.1007/s11214-024-01044-7)
  * [18] Pierel JDR, Arendse N, Ertl S, Huang X, Moustakas LA, Schuldt S, Shajib AJ, Shu Y, Birrer S, Bronikowski M, Hjorth J, Suyu SH, Agarwal S, Agnello A, Bolton AS, Chakrabarti S, Cold C, Courbin F, Della Costa JM, Dhawan S, Engesser M, Fox OD, Gall C, Gomez S, Goobar A, Jha SW, Jimenez C, Johansson J, Larison C, Li G, Marques-Chaves R, Mao S, Mazzali PA, Perez-Fournon I, Petrushevska T, Poidevin F, Rest A, Sheu W, Shirley R, Silver E, Storfer C, Strolger LG, Treu T, Wojtak R, Zenati Y. 2023 LensWatch. I. Resolved HST Observations and Constraints on the Strongly Lensed Type Ia Supernova 2022qmx (“SN Zwicky”). ApJ 948, 115. (10.3847/1538-4357/acc7a6)
  * [19] Mörtsell E, Johansson J, Dhawan S, Goobar A, Amanullah R, Goldstein DA. 2020 Lens modelling of the strongly lensed Type Ia supernova iPTF16geu. MNRAS 496, 3270–3280. (10.1093/mnras/staa1600)
  * [20] Diego JM, Bernstein G, Chen W, Goobar A, Johansson JP, Kelly PL, Mörtsell E, Nightingale JW. 2022 Microlensing and the type Ia supernova iPTF16geu. A&A 662, A34. (10.1051/0004-6361/202143009)
  * [21] Dhawan S, Johansson J, Goobar A, Amanullah R, Mörtsell E, Cenko SB, Cooray A, Fox O, Goldstein D, Kalender R, Kasliwal M, Kulkarni SR, Lee WH, Nayyeri H, Nugent P, Ofek E, Quimby R. 2020 Magnification, dust, and time-delay constraints from the first resolved strongly lensed Type Ia supernova iPTF16geu. MNRAS 491, 2639–2654. (10.1093/mnras/stz2965)
  * [22] Goldstein DA, Nugent PE, Goobar A. 2019 Rates and Properties of Supernovae Strongly Gravitationally Lensed by Elliptical Galaxies in Time-domain Imaging Surveys. ApJS 243, 6. (10.3847/1538-4365/ab1fe0)
  * [23] Sagués Carracedo A, Goobar A, Mörtsell E, Arendse N, Johansson J, Townsend A, Dhawan S, Nordin J, Sollerman J, Schulze S. 2024 Detectability and Characterisation of Strongly Lensed Supernova Lightcurves in the Zwicky Transient Facility. arXiv e-prints p. arXiv:2406.00052.
  * [24] Bellm EC, Kulkarni SR, Graham MJ, Dekany R, Smith RM, Riddle R, Masci FJ, Helou G, Prince TA, Adams SM, Barbarino C, Barlow T, Bauer J, Beck R, Belicki J, Biswas R, Blagorodnova N, Bodewits D, Bolin B, Brinnel V, Brooke T, Bue B, Bulla M, Burruss R, Cenko SB, Chang CK, Connolly A, Coughlin M, Cromer J, Cunningham V, De K, Delacroix A, Desai V, Duev DA, Eadie G, Farnham TL, Feeney M, Feindt U, Flynn D, Franckowiak A, Frederick S, Fremling C, Gal-Yam A, Gezari S, Giomi M, Goldstein DA, Golkhou VZ, Goobar A, Groom S, Hacopians E, Hale D, Henning J, Ho AYQ, Hover D, Howell J, Hung T, Huppenkothen D, Imel D, Ip WH, Ivezić Ž, Jackson E, Jones L, Juric M, Kasliwal MM, Kaspi S, Kaye S, Kelley MSP, Kowalski M, Kramer E, Kupfer T, Landry W, Laher RR, Lee CD, Lin HW, Lin ZY, Lunnan R, Giomi M, Mahabal A, Mao P, Miller AA, Monkewitz S, Murphy P, Ngeow CC, Nordin J, Nugent P, Ofek E, Patterson MT, Penprase B, Porter M, Rauch L, Rebbapragada U, Reiley D, Rigault M, Rodriguez H, van Roestel J, Rusholme B, van Santen J, Schulze S, Shupe DL, Singer LP, Soumagnac MT, Stein R, Surace J, Sollerman J, Szkody P, Taddia F, Terek S, Van Sistine A, van Velzen S, Vestrand WT, Walters R, Ward C, Ye QZ, Yu PC, Yan L, Zolkower J. 2019 The Zwicky Transient Facility: System Overview, Performance, and First Results. PASP 131, 018002. (10.1088/1538-3873/aaecbe)
  * [25] Fremling C, Miller AA, Sharma Y, Dugas A, Perley DA, Taggart K, Sollerman J, Goobar A, Graham ML, Neill JD, Nordin J, Rigault M, Walters R, Andreoni I, Bagdasaryan A, Belicki J, Cannella C, Bellm EC, Cenko SB, De K, Dekany R, Frederick S, Golkhou VZ, Graham MJ, Helou G, Ho AYQ, Kasliwal MM, Kupfer T, Laher RR, Mahabal A, Masci FJ, Riddle R, Rusholme B, Schulze S, Shupe DL, Smith RM, van Velzen S, Yan L, Yao Y, Zhuang Z, Kulkarni SR. 2020 The Zwicky Transient Facility Bright Transient Survey. I. Spectroscopic Classification and the Redshift Completeness of Local Galaxy Catalogs. ApJ 895, 32. (10.3847/1538-4357/ab8943)
  * [26] Sainz de Murieta A, Collett TE, Magee MR, Weisenbach L, Krawczyk CM, Enzi W. 2023 Lensed Type Ia supernovae in light of SN Zwicky and iPTF16geu. MNRAS 526, 4296–4307. (10.1093/mnras/stad3031)
  * [27] Townsend A, Nordin J, Sagués Carracedo A, Kowalski M, Arendse N, Dhawan S, Goobar A, Johansson J, Mörtsell E, Schulze S, Andreoni I, Fernández E, Kim AG, Nugent PE, Prada F, Rigault M, Sarin N, Sharma D, Bellm EC, Coughlin MW, Dekany R, Groom SL, Lacroix L, Laher RR, Riddle R, Aguilar J, Ahlen S, Bailey S, Brooks D, Claybaugh T, de la Macorra A, Dey A, Dey B, Doel P, Fanning K, Forero-Romero JE, Gaztañaga E, Gontcho SGA, Honscheid K, Howlett C, Kisner T, Kremin A, Lambert A, Landriau M, Le Guillou L, Levi ME, Manera M, Meisner A, Miquel R, Moustakas J, Mueller E, Myers AD, Nie J, Palanque-Delabrouille N, Poppett C, Rezaie M, Rossi G, Sanchez E, Schlegel D, Schubnell M, Seo H, Sprayberry D, Tarlé G, Zou H. 2024 Candidate strongly-lensed Type Ia supernovae in the Zwicky Transient Facility archive. arXiv e-prints p. arXiv:2405.18589. (10.48550/arXiv.2405.18589)
  * [28] Kamionkowski M, Riess AG. 2023 The Hubble Tension and Early Dark Energy. Annual Review of Nuclear and Particle Science 73, 153–180. (10.1146/annurev-nucl-111422-024107)
  * [29] Birrer S, Millon M, Sluse D, Shajib AJ, Courbin F, Koopmans LVE, Suyu SH, Treu T. 2022 Time-Delay Cosmography: Measuring the Hubble Constant and other cosmological parameters with strong gravitational lensing. arXiv e-prints p. arXiv:2210.10833. (10.48550/arXiv.2210.10833)
  * [30] Falco EE, Gorenstein MV, Shapiro II. 1985 On model-dependent bounds on H 0 from gravitational images : application to Q 0957+561 A, B.. ApJ 289, L1–L4. (10.1086/184422)
  * [31] Birrer S, Dhawan S, Shajib AJ. 2021 The Hubble constant from strongly lensed supernovae with standardizable magnifications. arXiv e-prints p. arXiv:2107.12385.
  * [32] Rodney SA, Strolger LG, Kelly PL, Bradač M, Brammer G, Filippenko AV, Foley RJ, Graur O, Hjorth J, Jha SW, McCully C, Molino A, Riess AG, Schmidt KB, Selsing J, Sharon K, Treu T, Weiner BJ, Zitrin A. 2016 SN Refsdal: Photometry and Time Delay Measurements of the First Einstein Cross Supernova. ApJ 820, 50. (10.3847/0004-637X/820/1/50)
  * [33] Kelly PL, Rodney S, Treu T, Oguri M, Chen W, Zitrin A, Birrer S, Bonvin V, Dessart L, Diego JM, Filippenko AV, Foley RJ, Gilman D, Hjorth J, Jauzac M, Mandel K, Millon M, Pierel J, Sharon K, Thorp S, Williams L, Broadhurst T, Dressler A, Graur O, Jha S, McCully C, Postman M, Schmidt KB, Tucker BE, von der Linden A. 2023 Constraints on the Hubble constant from supernova Refsdal’s reappearance. Science 380, abh1322. (10.1126/science.abh1322)
  * [34] Chen W, Kelly PL, Frye BL, Pierel J, Willner SP, Pascale M, Cohen SH, Conselice CJ, Engesser M, Furtak LJ, Gilman D, Grogin NA, Huber S, Jha SW, Johansson J, Koekemoer AM, Larison C, Meena AK, Siebert MR, Windhorst RA, Yan H, Zitrin A. 2024 JWST Spectroscopy of SN H0pe: Classification and Time Delays of a Triply-imaged Type Ia Supernova at z = 1.78. arXiv e-prints p. arXiv:2403.19029. (10.48550/arXiv.2403.19029)
  * [35] Pierel JDR, Frye BL, Pascale M, Caminha GB, Chen W, Dhawan S, Gilman D, Grayling M, Huber S, Kelly P, Thorp S, Arendse N, Birrer S, Bronikowski M, Canameras R, Coe D, Cohen SH, Conselice CJ, Driver SP, Dsilva JCJ, Engesser M, Foo N, Gall C, Garuda N, Grillo C, Grogin NA, Henderson J, Hjorth J, Jansen RA, Johansson J, Kamieneski PS, Koekemoer AM, Larison C, Marshall MA, Moustakas LA, Nonino M, Ortiz, R. I, Petrushevska T, Pirzkal N, Robotham A, Ryan, R. E. J, Schuldt S, Strolger LG, Summers J, Suyu SH, Treu T, Willmer CNA, Windhorst RA, Yan H, Zitrin A, Acebron A, Chakrabarti S, Coulter DA, Fox OD, Huang X, Jha SW, Li G, Mazzali PA, Meena AK, Perez-Fournon I, Poidevin F, Rest A, Riess AG. 2024 JWST Photometric Time-Delay and Magnification Measurements for the Triply-Imaged Type Ia “Supernova H0pe” at z = 1.78. arXiv e-prints p. arXiv:2403.18954. (10.48550/arXiv.2403.18954)
  * [36] Pascale M, Frye BL, Pierel JDR, Chen W, Kelly PL, Cohen SH, Windhorst RA, Riess AG, Kamieneski PS, Diego JM, Meena AK, Cha S, Oguri M, Zitrin A, Jee MJ, Foo N, Leimbach R, Koekemoer AM, Conselice CJ, Dai L, Goobar A, Siebert MR, Strolger L, Willner SP. 2024 SN H0pe: The First Measurement of $H_{0}$ from a Multiply-Imaged Type Ia Supernova, Discovered by JWST. arXiv e-prints p. arXiv:2403.18902. (10.48550/arXiv.2403.18902)
  * [37] Hsiao EY, Conley A, Howell DA, Sullivan M, Pritchet CJ, Carlberg RG, Nugent PE, Phillips MM. 2007 K-Corrections and Spectral Templates of Type Ia Supernovae. ApJ 663, 1187–1200. (10.1086/518232)
  * [38] Johansson J, Goobar A, Price SH, Sagués Carracedo A, Della Bruna L, Nugent PE, Dhawan S, Mörtsell E, Papadogiannakis S, Amanullah R, Goldstein D, Cenko SB, De K, Dugas A, Kasliwal MM, Kulkarni SR, Lunnan R. 2021 Spectroscopy of the first resolved strongly lensed Type Ia supernova iPTF16geu. MNRAS 502, 510–520. (10.1093/mnras/staa3829)
  * [39] Pierel JR, Rodney SA. 2019 SNTD: Supernova Time Delays. Astrophysics Source Code Library, record ascl:1902.001.
  * [40] Guy J, Astier P, Baumont S, Hardin D, Pain R, Regnault N, Basa S, Carlberg RG, Conley A, Fabbro S, Fouchez D, Hook IM, Howell DA, Perrett K, Pritchet CJ, Rich J, Sullivan M, Antilogus P, Aubourg E, Bazin G, Bronder J, Filiol M, Palanque-Delabrouille N, Ripoche P, Ruhlmann-Kleider V. 2007 SALT2: using distant supernovae to improve the use of type Ia supernovae as distance indicators. A&A 466, 11–21.
  * [41] Amanullah R, Johansson J, Goobar A, Ferretti R, Papadogiannakis S, Petrushevska T, Brown PJ, Cao Y, Contreras C, Dahle H, Elias-Rosa N, Fynbo JPU, Gorosabel J, Guaita L, Hangard L, Howell DA, Hsiao EY, Kankare E, Kasliwal M, Leloudas G, Lundqvist P, Mattila S, Nugent P, Phillips MM, Sandberg A, Stanishev V, Sullivan M, Taddia F, Östlin G, Asadi S, Herrero-Illana R, Jensen JJ, Karhunen K, Lazarevic S, Varenius E, Santos P, Sridhar SS, Wallström SHJ, Wiegert J. 2015 Diversity in extinction laws of Type Ia supernovae measured between 0.2 and 2 $\mu$m. MNRAS 453, 3300–3328. (10.1093/mnras/stv1505)
  * [42] More A, Suyu SH, Oguri M, More S, Lee CH. 2017 Interpreting the Strongly Lensed Supernova iPTF16geu: Time Delay Predictions, Microlensing, and Lensing Rates. ApJ 835, L25. (10.3847/2041-8213/835/2/L25)
  * [43] Dobler G, Keeton CR. 2006 Microlensing of Lensed Supernovae. ApJ 653, 1391–1399. (10.1086/508769)
  * [44] Foxley-Marrable M, Collett TE, Vernardos G, Goldstein DA, Bacon D. 2018 The impact of microlensing on the standardization of strongly lensed Type Ia supernovae. MNRAS 478, 5081–5090. (10.1093/mnras/sty1346)
  * [45] Weisenbach L, Collett T, Sainz de Murieta A, Krawczyk C, Vernardos G, Enzi W, Lundgren A. 2024 How to Break the Mass Sheet Degeneracy with the Lightcurves of Microlensed Type Ia Supernovae. arXiv e-prints p. arXiv:2403.03264. (10.48550/arXiv.2403.03264)
  * [46] Wojtak R, Hjorth J, Gall C. 2019 Magnified or multiply imaged? - Search strategies for gravitationally lensed supernovae in wide-field surveys. MNRAS 487, 3342–3355. (10.1093/mnras/stz1516)
  * [47] Arendse N, Dhawan S, Sagués Carracedo A, Peiris HV, Goobar A, Wojtak R, Alves C, Biswas R, Huber S, Birrer S, The LSST Dark Energy Science Collaboration. 2023 Detecting strongly-lensed type Ia supernovae with LSST. arXiv e-prints p. arXiv:2312.04621. (10.48550/arXiv.2312.04621)


# StAnD: A Dataset of Linear Static Analysis Problems

Luca Grementieri
Zuru Tech
<EMAIL_ADDRESS>
&Francesco Finelli
Zuru Tech
<EMAIL_ADDRESS>

###### Abstract

Static analysis of structures is a fundamental step for determining the
stability of structures. Both linear and non-linear static analyses consist of
the resolution of sparse linear systems obtained by the finite element method.
The development of fast and optimized solvers for sparse linear systems
appearing in structural engineering requires data to compare existing
approaches, tune algorithms or to evaluate new ideas. We introduce the Static
Analysis Dataset (StAnD) containing 303.000 static analysis problems obtained
applying realistic loads to simulated frame structures. Along with the
dataset, we publish a detailed benchmark comparison of the running time of
existing solvers both on CPU and GPU. We release the code used to generate the
dataset and benchmark existing solvers on Github. To the best of our
knowledge, this is the largest dataset for static analysis problems and it is
the first public dataset of sparse linear systems (containing both the matrix
and a realistic constant term).

##  1 Introduction

Structural analysis represents the set of theories and methods which allow to
idealize a structure through a mathematical model able to predict their
response in terms of displacements and internal stresses, when subject to
various types of external actions. The determination of response quantities
plays a crucial role for the design of safe and robust constructions in civil
engineering applications.

The structural behavior of a solid body is analytically described by partial
differential equations derived from the continuum mechanics principles of
equilibrium and compatibility, related through the material constitutive
relations. Unfortunately, analytical models of structures rarely have closed
form solutions in practical engineering problems. Consequently, the continuum
system is usually approximated with a discrete system with finite number of
degrees of freedom, while the set of governing differential equations is
converted and approximated with a system of algebraic equations, solvable with
a variety of numerical methods.

In this context, the Finite Element Method represents the most widespread
technique for numerical analysis of structures. The method produces a discrete
idealization of the structure which in static domain leads to a sparse
algebraic system

$K\mathbf{u}=\mathbf{f}$

where the variable term $\mathbf{u}$, the constant term $\mathbf{f}$, and the
coefficient matrix K respectively represent the displacement vector, the
external load vector, and the global stiffness matrix.

The solution of such algebraic system provides the values of generalized
displacements at each degree of freedom of the structure consequent to the
prescribed loading conditions. Once displacements are known, then strains and
stresses can be computed following kinematic compatibility conditions and
material constitutive law, so that all response quantities necessary for
structural design procedures are determined.

Many algorithms for solving sparse linear systems are published at a great
pace, but currently there is no standard dataset to compare their running time
on real problems.

We introduce Static Analysis Dataset (StAnD) to formalize the evaluation of
new resolution methods and to spur research in resolution methods specifically
tailored to structural engineering and static analysis problems.

To best of our knowledge there is no existing large dataset of sparse linear
systems. A few sparse matrices (often less than 10) are published for many
engineering problems in the Matrix Market (Boisvert et al., 1997) or in the
SuiteSparse matrix collection (Davis and Hu, 2011), but their limited number
is not sufficient to measure the running time in the average case or to
measure reliably how the resolution algorithm scales with the size and the
number of non-zeros in the sparse matrix. Additionally, no constant term is
provided in conjunction with the matrices. The constant term, is maybe not
fundamental when direct methods are used, but it becomes important if we want
to measure the effectiveness of iterative methods, whose behavior depends on
the relationship between the initialization of the solution and the real
solution.

A large-scale dataset can also pave the way to novel methods based on machine
learning. Large-scale datasets have led to incredible advances in many areas
of artificial intelligence in the last few years, like computer vision and
natural language processing and we hope to foster research in machine learning
application for structural engineering releasing this dataset. For this reason
StAnD is already divided into a training split and a test split.

In summary, our contribution is threefold:

  * •

We publish a novel dataset of static analysis problems of frame structures.
The dataset is composed of a training set of 300.000 problems and 3.000 test
problems. Every split is evenly divided into
small,111https://storage.googleapis.com/zurutech-public-
datasets/stand/stand_small.zip
medium222https://storage.googleapis.com/zurutech-public-
datasets/stand/stand_medium.zip and large
problems.333https://storage.googleapis.com/zurutech-public-
datasets/stand/stand_large.zip This clear separation allow to inspect the
scaling behavior of solvers. The structures are grid-like structures
procedurally designed, while the loads applied are obtained simulating
elementary actions of permanent (proper weights of the structure and of non-
structural elements) and variable loads (furniture, people, wind, snow).

  * •

We conduct a thorough evaluation of existing open-source implementations of
direct and iterative methods for the resolution of sparse linear systems. In
particular, we put relevant attention into implementations that support GPUs
since GPUs are quite common nowadays, but many algorithms are not really
suited to leverage their processing power at its best.

  * •

We release the code to generate the dataset based on OpenSeesPy (Zhu et al.,
2018) and we publish the code and the Docker container used to benchmark
existing solvers.444https://github.com/zurutech/stand

##  2 Related Work

###  2.1 Existing Datasets

One the earliest collection of sparse matrices is the Matrix Market (Boisvert
et al., 1997), which includes in particular the Harwell-Boeing Sparse Matrix
Collection (Duff et al., 1989) and SPARSEKIT collection (Saad, 1994). At the
moment, the largest and most comprehensive collection of sparse matrices
publicly available is the SuiteSparse Matrix collection (Davis and Hu, 2011),
which includes also a large subset of the Matrix Market.

These collections are generic with sparse matrices from several domains like
aerodynamics, fluid dynamics, structural engineering and economics.
Additionally they are very small: in fact the SuiteSparse Matrix collection is
composed of less than 3000 matrices. For every single application, they
contain only a handful of matrices and it is difficult to understand if the
available data faithfully reproduces the setting of interest because no
metadata about the problem that has led to such matrices exists. These
datasets can be used to compare generic solvers across several domains, but
they are not so useful to compare algorithms for a specific domain of
application. Furthermore, both collections only contain matrices and if we
want to test an algorithm for the resolution of linear system, we have to
sample a random vector of constant terms. The sampled constant terms can be
physically unrealistic or infeasible and their choice can influence the
evaluation of iterative methods.

###  2.2 Methods

The number of methods for the resolution of a sparse linear system
$A\mathbf{x}=\mathbf{b}$ is large and many approaches have been suggested to
tackle this classical problem. Linear systems for the static analysis of
structures are always symmetric positive definite, thus we focus on algorithms
applicable in this setting. The methods can be classified into three groups:
direct, iterative and hybrid methods. We describe here the methods implemented
in the most used open source libraries because these are the methods included
in our benchmark.

####  2.2.1 Direct Methods

Davis et al. (2016) give a comprehensive review of direct methods and their
historical development. Most direct method work in three stages: matrix
permutation, matrix factorization, and triangular system resolution.

Matrix permutation. A sparse matrix must typically be permuted either before
or during its numeric factorization, either for reducing fill-in or for
numerical stability. _Fill-in_ is the introduction of new non-zeros in the
factors that do not appear in the corresponding positions in the matrix being
factorized. Finding the permutation that minimizes the number of fill-in is a
NP-hard problem (Rose and Tarjan, 1978; Yannakakis, 1981; Luce and Ng, 2014),
so the literature propose many heuristic fill-reducing methods. The libraries
included in our benchmark implement the following permutation algorithms:
symmetric reverse Cuthill-McKee ordering (RCM) (Cuthill and McKee, 1969),
minimum degree ordering (MMD) (George and Liu, 1989), symmetric approximate
minimum degree ordering (AMD) (Amestoy et al., 1996), column approximate
minimum degree ordering (COLAMD) (Davis et al., 2004), METIS (Karypis and
Kumar, 1998) and its parallel variant ParMETIS (LaSalle and Karypis, 2016).

Matrix factorization. There are three main families of factorizations used for
sparse linear system resolution: QR factorization ($A=QR$, with $Q$ an
orthogonal matrix and $R$ an upper triangular matrix), LU factorization
($A=LU$, with $L$ a lower triangular matrix and $U$ an upper triangular
matrix) and, if $A$ is symmetric positive-definite, Cholesky factorization
($A=LL^{T}$, with $L$ a lower triangular matrix). Several specialized
algorithms exist to perform such decompositions efficiently with sparse
matrices. In almost all methods, the factorization splits into two phases: a
symbolic phase that typically depends only on the non-zero pattern of $A$, and
a numerical phase that produces the factorization itself. The symbolic phase
finds the non-zero pattern of the triangular factors of $A$ without
determining the values of the triangular factors itself. This phase improves
the efficiency of the subsequent numerical phase and its result can be reused
on matrices with identical non-zero pattern, a common situation when solving
non-linear or differential equations. The numerical phase can be implemented
in many different variants, even parallel ones. We will refer to the specific
implementations in section 4.1.

Triangular system resolution. After the factorization of the sparse matrix
$A$, we can recover the solution solving one or two consecutive triangular
systems. Using the QR factorization, the triangular system to be solved is

$R\mathbf{x}=\mathbf{y}\quad\text{where}\quad y=Q^{T}\mathbf{b},$

while for the LU factorization (and similarly for Cholesky factorization) the
final solution is obtained solving two consecutive triangular systems

$L\mathbf{y}=\mathbf{b}\quad\text{then}\quad U\mathbf{x}=\mathbf{y}.$

The algorithm to solve (sparse) triangular systems is quite similar to the
matrix-vector multiplication algorithm, but it is inherently sequential. This
implies that it is efficient like matrix-vector multiplication on CPUs, but on
GPUs it is much slower.

####  2.2.2 Iterative Methods

Iterative methods are increasingly popular for the resolution of sparse linear
systems. Here we list only methods implemented in the open-source libraries of
our benchmark. Please refer to the book by Saad (2003) for a detailed and
complete introduction. Iterative methods are the workhorse of modern hybrid
methods, so we are going to introduce them together. Hybrid methods apply an
existing iterative algorithm to a preconditioned sparse system. The
preconditioning transforms the linear system in an equivalent one, while
reducing the condition number of the sparse matrix. Since the condition number
determines the rate of convergence of iterative algorithms, this preliminary
step can greatly reduce the number of iterations needed and consequently the
running time of the algorithm.

Krylov methods. The most successful iterative methods are the conjugate
gradient (CG) for symmetric positive-definite matrices and the generalized
minimal residual method (GMRES) (Saad and Schultz, 1986) for non-symmetric
matrices. Both methods fall into the family of Krylov methods since they are
based on projections onto Krylov subspaces. A _Krylov subspace_ is the
subspace spanned by vectors of the form $p(A)b$, where $p$ is a polynomial. In
short, these techniques approximate $A^{-1}\mathbf{b}$ by $p(A)\mathbf{b}$.
Other methods included in our benchmark are flexible conjugate gradient
(Notay, 2000), flexible GMRES (Giraud et al., 2010), quasi minimal residual
(Chan et al., 1998), conjugate residual method (Hestenes and Stiefel, 1952),
minimal residual and symmetric LQ (Paige and Saunders, 1975).

Preconditioning. A _preconditioner_ $P$ of a matrix $A$ is a matrix such that
$P^{-1}A$ or $AP^{-1}$ have a smaller condition number than $A$. The _left
preconditioned_ system is the equivalent system

$P^{-1}A\mathbf{x}=P^{-1}\mathbf{b}.$

Instead, _right preconditioning_ requires the resolution of two systems:

$AP^{-1}\mathbf{y}=\mathbf{b}\text{ then }P\mathbf{x}=\mathbf{y}.$

Most preconditioners are inspired from direct methods or classic iterative
methods. For example, the incomplete factorization ILU(0) decomposes the
matrix $A$ as $A=LU-R$ where $L$ and $U$ have the same nonzero structure as
the lower and upper parts of $A$ respectively, and $R$ is the residual or
error of the factorization. The preconditioners covered by the libraries in
our benchmark are Jacobi, SOR, SSOR, ILU(0), BlockILU(0), ICC.

##  3 Dataset Description

The StAnD dataset is composed of 303.000 static analysis problems of frame
structures divided into 6 parts: 100.000 small training problems, 100.000
medium training problems, 100.000 large training problems, 1.000 small test
problems, 1.000 medium test problems and 1.000 large training problems. The
size of a problem is determined by the number of degrees of freedom (DOFs) of
the structure model or equivalently by the number of rows or columns of the
stiffness matrix associated to the structure. Small problems have 2115 DOFs on
average (and at most 5166 DOFs), medium problems have about 7000 DOFs (and at
most 14718 DOFs), while large problems can have up to 31770 DOFs with about
15500 DOFs on average. The division of problems by size allows to compare the
scalability of solvers.

The dataset is programmatically created using OpenSeesPy (Zhu et al., 2018),
the Python interface to the OpenSees finite element solver (Mazzoni et al.,
2006). For each structural model with its loading configuration, a static
analysis is performed in order to compute nodal displacements. Therefore,
every problem in the dataset is a tuple $(K,\mathbf{f},\mathbf{u})$, where $K$
is the sparse stiffness matrix associated to the structure, $\mathbf{f}$ is a
load vector and $\mathbf{u}$ is the ground-truth displacement vector such that
$K\mathbf{u}=\mathbf{f}$. In the training set, for the same structure (i.e.
the same stiffness matrix $K$), we apply different load configurations. In the
test set, we use a single load configuration for every structure to maximize
the variability.

###  3.1 Frame Structure Generation

(a)

(b)

(c)

Figure 1: A frame structure is created from a regular 3D grid (a), removing
cubes (b) and then changing the length of beams (c).

Structures are created from a regular 3D grid composed of $B\times D\times H$
unit cubes. The values $B,D,H$ that determine the size of the starting grid
are sampled uniformly in ranges according to the category (small, medium, or
large) of the structure. From the regular grid, we remove a random number of
cubes (from 0 to $H$ included) from each of the $b\times d$ columns made of
$h$ cubes. The cubes are removed from the top, creating structures composed of
stacks of cubes of different heights. Then the carved grid and its cubes are
deformed: the thickness of every plane of cubes is changed from 1m to a random
value between 3 and 6 meters. This applies to horizontal planes of cubes
orthogonal to the $z$ axis (i.e. with cubes at the same level of height) and
to vertical planes orthogonal to the $x$ or the $y$ axis. This deformation
transforms the cubes into rectangular parallelepipeds. Finally, the frame
structure is created from the edges of the deformed grid. The deformation of
the grid is instrumental to vary the lengths of beams.

The beams elements used follow the Timoshenko beam model. The parameters
needed to fully characterize the beam model are sampled from the ranges shown
in Table 1 and they are shared by all beams in the same structure.

Table 1: Beam parameters ranges Length (m) | Width (m) | Height (m) | Young’s modulus (Pa) | Density (kg/$\text{m}^{3}$)
---|---|---|---|---
$[3.0,6.0]$ | $[0.2,0.4]$ | $[0.4,0.6]$ | $[2.5\cdot 10^{10},4.3\cdot 10^{10}]$ | $[2.0\cdot 10^{3},2.6\cdot 10^{3}]$

###  3.2 Load Application

The verification of a structure consists in its verification against ultimate
limit state (ULS) and serviceability limit state (SLS). The loading conditions
at limit states are obtained defining various combinations of several
elementary loads. The elementary actions are the ones induced by various
sources of gravity loads, such as permanent structural and non-structural dead
weights or variable live loads, by the atmospheric actions, such as wind
pressure, or by earthquakes. In linear static analysis, we do not need to
solve a sparse linear system for every combination of such actions, because
the principle of superposition allows to solve the algebraic system associated
to each elementary condition separately and then linearly combine the
solutions. For example, if a structure’s ULS is obtained as the combination of
the dead weight of the structure, a variable live load and the action of wind,
then the load vector is

$\mathbf{f}=\varphi_{g}\mathbf{f}_{g}+\varphi_{v}\mathbf{f}_{v}+\varphi_{w}\mathbf{f}_{w}.$

In this case, the displacement vector $\mathbf{u}$ such that
$K\mathbf{u}=\mathbf{f}$ is

$\mathbf{u}=\varphi_{g}\mathbf{u}_{g}+\varphi_{v}\mathbf{u}_{v}+\varphi_{w}\mathbf{u}_{w}\text{
where
}K\mathbf{u}_{g}=\mathbf{f}_{g},K\mathbf{u}_{v}=\mathbf{f}_{v},K\mathbf{u}_{w}=\mathbf{f}_{w}.$

Following the engineering practice, in StAnD we include only problems obtained
applying a single elementary action. The considered load typologies have been
evaluated and computed according to the prescriptions of the Italian building
code (Ministero delle Infrastrutture e dei Trasporti, 2018).

Proper weight of the structure. For every structural member in StAnD we sample
a density $\rho$, a section width $S_{w}$ and a section height $S_{h}$ from
the ranges specified in Table 1. We denote with $g$ the standard gravitational
acceleration for the surface of the Earth. For beams the proper weight load is
a uniform load directed towards the negative direction of the $z$ axis with
magnitude

$p_{g}=g\cdot\rho\cdot S_{w}\cdot S_{h}.$

Variable vertical pressure. In addition to its proper weight, the structure
has to support the weight of slabs and additional variable loads linked to the
function of the structure itself. Slabs could be modelled more precisely with
a 4-node element, but since we are working with frame structures only, we have
to compute the nodal equivalent loads. We assume that every slab transfers all
loads to the horizontal beams just below it. Under this simplifying
assumption, a uniform load on a rectangular slab is equivalent to a pair of
identical triangle loads on the shortest sides of the rectangle and a pair of
trapezoidal loads on the longest sides of the rectangle, as shown in Figure In
practice, triangle and trapezoidal loads are further simplified into uniform
loads with the same total force on the beam. In conclusion, denoting $\ell$
and $L$ the short and long dimensions of the beams of the slab, the uniform
loads on beams of the slab have magnitude

$p_{v}^{\text{short}}=P\cdot\frac{\ell}{4}\qquad
p_{v}^{\text{long}}=P\cdot\frac{\ell}{2L}\cdot\left(L-\frac{\ell}{2}\right),$

where $P$ is the vertical pressure on the slab.

Wind pressure. The action of the wind is modelled as pressures acting on the
windward face of the structure and depressions on the leeward and lateral
faces. The direction of the wind in static analysis is taken parallel to the
principal axes of the building plan. The pressure of the wind depends on the
altitude $z$ and it is given by

$p_{w}(z)=\frac{1}{2}\rho_{\text{air}}\cdot v^{2}\cdot c_{e}(z)\cdot c_{p},$

where $\rho_{air}=1.25\text{kg/m}^{3}$ is the density of the air, $v$ is the
speed of wind, $c_{e}(z)$ is the coefficient of exposition (which is a non-
decreasing function of $z$) and $c_{p}$ is the pressure coefficient of the
face under examination. The coefficient of exposition is defined as

$c_{e}(z)=\begin{cases}\displaystyle
k^{2}\ln\left(\frac{z}{z_{0}}\right)\left[7+\ln\left(\frac{z}{z_{0}}\right)\right]&\text{
if }z\geq z_{\text{min}}\\\ c_{e}(z_{min})&\text{ if
}z<z_{\text{min}}.\end{cases}$

The parameters $k$, $z_{0}$ and $z_{\text{min}}$ depend on the category of
exposition of the building site and they are chosen according to Table 2.

Table 2: Parameters determining the coefficient of exposition Category of exposition | $k$ |  $z_{0}$ (m) |  $z_{\text{min}}$ (m)
---|---|---|---
I | $0.17$ | $0.01$ | $2$
II | $0.19$ | $0.05$ | $4$
III | $0.20$ | $0.10$ | $5$
IV | $0.22$ | $0.30$ | $8$
V | $0.23$ | $0.70$ | $12$

The pressure coefficient depends on the ratio between the height of the
structure $h$ and its depth $d$ measured along the axis parallel to the
direction of the wind. The positive sign of the pressure coefficient indicates
a pressure, while the negative sign signal a depression. The pressure
coefficient has different formulations for the windward, the leeward and the
lateral faces of the structure:

$c_{p}^{\text{windward}}(h,d)=\begin{cases}0.7+0.1\frac{h}{d}&\text{ if
}\frac{h}{d}\leq 1\\\ 0.8&\text{ if }\frac{h}{d}>1;\end{cases}$
$c_{p}^{\text{leeward}}(h,d)=\begin{cases}-0.3-0.2\frac{h}{d}&\text{ if
}\frac{h}{d}\leq 1\\\ -0.5-0.05(\frac{h}{d}-1)&\text{ if
}\frac{h}{d}>1;\end{cases}$
$c_{p}^{\text{lateral}}(h,d)=\begin{cases}-0.5-0.8\frac{h}{d}&\text{ if
}\frac{h}{d}\leq\frac{1}{2}\\\ -0.9&\text{ if
}\frac{h}{d}>\frac{1}{2}.\end{cases}$

We suppose that every external frame of structure (composed of a pair of
horizontal beams and a pair of vertical columns) is covered by a wall able to
transfer all the pressure of the wind on beams and columns. The load on a wall
is not uniform since it varies with the altitude, but for simplicity we fix
the value of the coefficient of exposition to $c_{e}(z^{\text{top}})$ fro the
entire loaded surface obtaining a uniform load. Since $c_{e}(z)$ is a non-
decreasing function and $z^{\text{top}}$ is the highest altitude of the wall,
the resulting load is an upper bound over the real load. Lastly, the wind
uniform load is distributed on beams using the same approximation introduced
for vertical pressure on slabs.

Snow pressure. The pressure of the snow on the roof of the structure is
modelled similarly to a variable vertical pressure, but it can only be applied
to slabs at the top of the structure. The formula modelling the snow pressure
is

$p_{s}=p_{0}\cdot\mu,$

where $p_{0}$ is the load of the snow on the ground and $\mu$ is a coefficient
that depends on the slope of the roof. The value $p_{0}$ is sampled in the
interval $[600.0,5600.0]$ to cover different weather conditions. The
coefficient $\mu=0.8$ is constant for all structures in the dataset since the
roof is always flat.

##  4 Benchmark

We use StAnD to measure the running time of linear solvers on static analysis
problem. Our benchmark try to be one of the most comprehensive, with a
particular focus on libraries with GPU support. We will see that not all
algorithms are faster on GPUs than on CPUs as we would expect. GPUs shines for
computations with a high degree of parallelization and low degree of inter-
communication. If the algorithm is inherently sequential and it does not
satisfy thee requirements, then the CPU’s running time can be lower than the
GPU’s running time.

###  4.1 Open-source Libraries

Our benchmark includes the following open-source C++ libraries with GPU
support: SuiteSparse (Davis, 2019), cuSOLVER (NVIDIA Corporation, 2021),
SuperLU_DIST (Li and Demmel, 2003), PETSc (Balay et al., 2021) and ViennaCL
(Rupp et al., 2016). SuiteSparse, cuSOLVER and SuperLU_DIST implements direct
methods, while PETSc and ViennaCL are focused on iterative methods. Since
PETSc and SuiteSparse do not fully support GPUs and they have better
performances on CPUs, we report only their results obtained on CPUs.

###  4.2 Results

We test the libraries on a PC equipped with 12 CPU cores (Intel i7-6850K CPU @
3.60GHz) and a GPU Nvidia GeForce GTX 1080 Ti. We measure the time needed to
solve the linear system, without considering the time spent to convert the
system in the format accepted by each library or the time needed to transfer
data on the GPU. For iterative algorithms we stop the iteration when the
relative error goes below $10^{-6}$ or when we reach 1000 iterations. Since we
known the OpenSees solution $x$, we can also evaluate the _standard error_ of
the approximation $\widetilde{x}$ with respect to $x$ defined as

$SE(\widetilde{x},x)=\frac{||x-\widetilde{x}||}{||x||}.$

We consider that an algorithm has converged for all problems if the average
standard error is below $10^{-3}$. When an algorithm does not converge, its
running time is not reported.

Tables 3 and 4 show the average running time of solvers for the test split of
StAnD on small, medium and large problems. As expected, the methods tailored
to solve SPD systems (Cholesky factorization, conjugate gradient and conjugate
residuals) are faster than methods specifically studied to solve non-SPD
systems (LU and QR factorizations, GMRES, QMR, Minimum Residual and SymmLQ).
We can notice that iterative methods are much faster than concurrent methods,
and preconditioning definitely helps convergence. In particular incomplete
factorization methods like ICC and ILU(0) gives a relevant speed-up on CPU
(from 4x to 25x). On GPU, preconditioning reduces the average standard error
but it often results in slower running times, probably because the
factorization cannot be easily parallelized. The same difficulty in
parallelization explains why direct algorithms are always faster on CPU than
on GPU.

In conclusion, we advice practitioners to use the (preconditioned) conjugate
gradient algorithm to solve linear static analysis problems. Comparing the
conjugate gradient algorithm implementations on CPU (PETSc) and on GPU
(ViennaCL), we see that the GPU is noticeably faster only on large problems:
in fact, for small and medium problems the GPU is severely underutilized. For
this reason, we recommend to use PETSc implementation with ILU(0)
preconditioning for analyses with less than 7000 DOFs and to optionally switch
to ViennaCL implementation on GPU for larger structures.

Table 3: Benchmark timing results on CPU. Library | Method |  Permutation / Preconditioner |  Small problems time (ms) |  Medium problems time (ms) |  Large problems time (ms)
---|---|---|---|---|---
SuiteSparse (CPU) | KLU | - | 44.8 | 600 | 3925
LU | - | 26.2 | 127 | 420
Cholesky | Simplicial | 21.9 | 227 | 896
Supernodal | 13.4 | 64.6 | 239
PETSc (CPU) |  Conjugate Gradient | - | 11.7 | 73.4 | 253
Jacobi | 9.35 | 63.1 | 229
SOR | 8.15 | 51.9 | 180
SSOR | 5.87 | 38.1 | 124
ICC | 2.28 | 17.0 | 53.1
ILU(0) | 2.13 | 15.8 | 48.6
Conjugate Residuals | - | 11.0 | 67.6 | 239
Jacobi | 8.30 | 53.9 | 194
SOR | 7.55 | 46.7 | 159
SSOR | 5.70 | 35.6 | 113
ICC | 2.33 | 16.8 | 51.8
ILU(0) | 2.18 | 15.7 | 47.9
GMRES | - | 74.7 | 149 | 442
Jacobi | 32.1 | 141 | 440
SOR | 16.7 | 105 | 407
SSOR | 13.9 | 85.1 | 305
ICC | 3.44 | 23.1 | 75.5
ILU(0) | 3.25 | 21.9 | 71.0
Flexible GMRES | - | 74.2 | 149 | 443
Jacobi | 38.3 | 146 | 447
SOR | 16.8 | 105 | 407
SSOR | 12.2 | 76.7 | 280
ICC | 3.13 | 21.2 | 69.8
ILU(0) | 2.96 | 20.2 | 66.0
Minimum Residual | - | 14.4 | 88.2 | 312
Jacobi | 11.7 | 74.9 | 271
SOR | 8.95 | 55.3 | 191
SSOR | 6.81 | 42.2 | 138
ICC | 2.53 | 18.1 | 57.1
ILU(0) | 2.38 | 17.1 | 53.3
QMR | - | 142 | 254 | 876
Jacobi | - | - | -
SOR | 63.6 | 201 | 740
SSOR | 54.7 | 155 | 529
ICC | 12.0 | 62.6 | 233
ILU(0) | 11.9 | 57.9 | 212
SymmLQ | - | 16.0 | 99.4 | 338
Jacobi | - | - | -
SOR | - | - | -
SSOR | - | - | -
ICC | - | - | -
ILU(0) | - | - | -
Cholesky | - | 63.7 | 721 | 3526
LU | - | 44.3 | 490 | 2811
QR | - | 65.5 | 331 | 1110

Table 4: Benchmark timing results on GPU. Library | Method |  Permutation / Preconditioner |  Small problems time (ms) |  Medium problems time (ms) |  Large problems time (ms)
---|---|---|---|---|---
cuSOLVER (GPU) | Cholesky | - | 40.4 | 310 | 1339
METIS | 30.7 | 136 | 393
SymAMD | 35.1 | 212 | 817
SymRCM | 34.5 | 226 | 866
QR | - | 103 | 806 | 4151
METIS | 154 | 1521 | 7466
SymAMD | 212 | 2974 | 20784
SymRCM | 89.1 | 594 | 2571
SuperLU (GPU) | LU | - | 1767 | 21399 | -
ColAMD | 418 | 775 | 1631
METIS | 374 | 531 | 835
Minimum Degree | 425 | 626 | 906
ViennaCL (GPU) |  Conjugate Gradient | - | 10.3 | 15.0 | 23.8
Row Scaling | 19.3 | 30.3 | 41.7
Jacobi | 19.1 | 29.9 | 41.0
Chow-Patel | 26.0 | 43.9 | 70.2
ICC | 6.26 | 19.0 | 43.5
ILU(0) | 6.58 | 19.4 | 43.6
Block ILU(0) | - | - | -
GMRES | - | - | - | -
Row Scaling | - | - | -
Jacobi | - | - | -
Chow-Patel | - | - | -
ICC | - | - | -
ILU(0) | - | - | -
Block ILU(0) | - | - | -

##  5 Conclusion

We introduce the Static Analysis Dataset (StAnD), a novel dataset of linear
static analysis problems. The dataset gives the possibility to evaluate direct
and iterative methods for the resolution of sparse symmetric positive-definite
linear systems. Such linear systems are obtained simulating complex frame
structures of various sizes and applying realistic loads to them.

Many open source libraries with partial or complete GPU support implementing
direct or iterative methods were thoroughly evaluated on this dataset. The
evaluations provide a first benchmark on this dataset and show that there is
still considerable room for improvement to leverage the parallel processing
power of GPUs, since in many settings CPUs are still faster than GPUs. It is
our hope that the proposed dataset will stimulate the development of new
sparse linear solvers specifically tailored for GPUs and for structural
analysis.

## References

  * Amestoy et al. [1996] Patrick R Amestoy, Timothy A Davis, and Iain S Duff.  An approximate minimum degree ordering algorithm.  _SIAM Journal on Matrix Analysis and Applications_ , 17(4):886–905, 1996.
  * Balay et al. [2021] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang.  PETSc/TAO users manual.  Technical Report ANL-21/39 - Revision 3.16, Argonne National Laboratory, 2021.
  * Boisvert et al. [1997] Ronald F Boisvert, Roldan Pozo, Karin Remington, Richard F Barrett, and Jack J Dongarra.  Matrix market: a web resource for test matrix collections.  In _Quality of Numerical Software_ , pages 125–137. Springer, 1997\.
  * Chan et al. [1998] Tony F Chan, Lisette De Pillis, and Henk Van Der Vorst.  Transpose-free formulations of lanczos-type methods for nonsymmetric linear systems.  _Numerical Algorithms_ , 17(1):51–66, 1998.
  * Cuthill and McKee [1969] Elizabeth Cuthill and James McKee.  Reducing the bandwidth of sparse symmetric matrices.  In _Proceedings of the 1969 24th national conference_ , pages 157–172, 1969.
  * Davis [2019] Timothy A Davis.  Algorithm 1000: Suitesparse: Graphblas: Graph algorithms in the language of sparse linear algebra.  _ACM Transactions on Mathematical Software (TOMS)_ , 45(4):1–25, 2019.
  * Davis and Hu [2011] Timothy A. Davis and Yifan Hu.  The university of florida sparse matrix collection.  _ACM Trans. Math. Softw._ , 38(1), 2011.  ISSN 0098-3500.  URL https://doi.org/10.1145/2049662.2049663.
  * Davis et al. [2004] Timothy A Davis, John R Gilbert, Stefan I Larimore, and Esmond G Ng.  A column approximate minimum degree ordering algorithm.  _ACM Transactions on Mathematical Software (TOMS)_ , 30(3):353–376, 2004.
  * Davis et al. [2016] Timothy A. Davis, Sivasankaran Rajamanickam, and Wissam M. Sid-Lakhdar.  A survey of direct methods for sparse linear systems.  _Acta Numerica_ , 25:383–566, 2016.
  * Duff et al. [1989] Iain S Duff, Roger G Grimes, and John G Lewis.  Sparse matrix test problems.  _ACM Transactions on Mathematical Software (TOMS)_ , 15(1):1–14, 1989.
  * George and Liu [1989] Alan George and Joseph WH Liu.  The evolution of the minimum degree ordering algorithm.  _Siam review_ , 31(1):1–19, 1989.
  * Giraud et al. [2010] Luc Giraud, Serge Gratton, Xavier Pinel, and Xavier Vasseur.  Flexible gmres with deflated restarting.  _SIAM Journal on Scientific Computing_ , 32(4):1858–1878, 2010.
  * Hestenes and Stiefel [1952] Magnus R Hestenes and Eduard Stiefel.  Methods of conjugate gradients for solving linear systems’.  _Journal of Research of the National Bureau of Standards_ , 49(6), 1952.
  * Karypis and Kumar [1998] George Karypis and Vipin Kumar.  A fast and high quality multilevel scheme for partitioning irregular graphs.  _SIAM Journal on scientific Computing_ , 20(1):359–392, 1998.
  * LaSalle and Karypis [2016] Dominique LaSalle and George Karypis.  A parallel hill-climbing refinement algorithm for graph partitioning.  In _2016 45th International Conference on Parallel Processing (ICPP)_ , pages 236–241. IEEE, 2016.
  * Li and Demmel [2003] Xiaoye S. Li and James W. Demmel.  SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems.  _ACM Trans. Mathematical Software_ , 29(2):110–140, June 2003.
  * Luce and Ng [2014] Robert Luce and Esmond G Ng.  On the minimum flops problem in the sparse cholesky factorization.  _SIAM Journal on Matrix Analysis and Applications_ , 35(1):1–21, 2014.
  * Mazzoni et al. [2006] Silvia Mazzoni, Frank McKenna, Michael H Scott, Gregory L Fenves, et al.  Opensees command language manual.  _Pacific Earthquake Engineering Research (PEER) Center_ , 264:137–158, 2006.
  * Ministero delle Infrastrutture e dei Trasporti [2018] Ministero delle Infrastrutture e dei Trasporti.  D.M. 17/01/2018. Aggiornamento delle «Norme tecniche perle costruzioni».  Gazzetta Ufficiale Serie Generale n.42 del 20/02/2018 - Suppl. Ordinario n. 8., 2018.
  * Notay [2000] Yvan Notay.  Flexible conjugate gradients.  _SIAM Journal on Scientific Computing_ , 22(4):1444–1460, 2000.
  * NVIDIA Corporation [2021] NVIDIA Corporation.  _cuSOLVER_.  NVIDIA Corporation, 2021.  URL https://docs.nvidia.com/cuda/cusolver/.
  * Paige and Saunders [1975] Christopher C Paige and Michael A Saunders.  Solution of sparse indefinite systems of linear equations.  _SIAM journal on numerical analysis_ , 12(4):617–629, 1975.
  * Rose and Tarjan [1978] Donald J Rose and Robert Endre Tarjan.  Algorithmic aspects of vertex elimination on directed graphs.  _SIAM Journal on Applied Mathematics_ , 34(1):176–197, 1978.
  * Rupp et al. [2016] Karl Rupp, Philippe Tillet, Florian Rudolf, Josef Weinbub, Andreas Morhammer, Tibor Grasser, Ansgar Jungel, and Siegfried Selberherr.  Viennacl—linear algebra library for multi-and many-core architectures.  _SIAM Journal on Scientific Computing_ , 38(5):S412–S439, 2016.
  * Saad and Schultz [1986] Youcef Saad and Martin H Schultz.  Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems.  _SIAM Journal on scientific and statistical computing_ , 7(3):856–869, 1986.
  * Saad [1994] Yousef Saad.  Sparskit: a basic tool kit for sparse matrix computations, 1994.
  * Saad [2003] Yousef Saad.  _Iterative methods for sparse linear systems_.  SIAM, 2003.
  * Yannakakis [1981] Mihalis Yannakakis.  Computing the minimum fill-in is np-complete.  _SIAM Journal on Algebraic Discrete Methods_ , 2(1):77–79, 1981.
  * Zhu et al. [2018] Minjie Zhu, Frank McKenna, and Michael H Scott.  Openseespy: Python library for the opensees finite element framework.  _SoftwareX_ , 7:6–11, 2018.