The LAPW method with eigendecomposition based on the Hari--Zimmermann generalized hyperbolic SVD

TL;DR

This paper introduces a highly parallel and accurate algorithm for the generalized eigendecomposition of complex Hermitian matrix pairs, leveraging a hyperbolic SVD approach within the context of Density Functional Theory.

Contribution

It presents a novel eigendecomposition algorithm based on the Hari--Zimmermann hyperbolic SVD, optimized for parallel computation and applicable to quantum mechanical Hamiltonians.

Findings

Algorithm achieves high accuracy in eigendecomposition

Parallel implementation improves computational efficiency

Applicable to complex Hermitian matrices in quantum physics

Abstract

In this paper we propose an accurate, highly parallel algorithm for the generalized eigendecomposition of a matrix pair $(H, S)$, given in a factored form $(F^{\ast} J F, G^{\ast} G)$. Matrices $H$ and $S$ are generally complex and Hermitian, and $S$ is positive definite. This type of matrices emerges from the representation of the Hamiltonian of a quantum mechanical system in terms of an overcomplete set of basis functions. This expansion is part of a class of models within the broad field of Density Functional Theory, which is considered the golden standard in condensed matter physics. The overall algorithm consists of four phases, the second and the fourth being optional, where the two last phases are computation of the generalized hyperbolic SVD of a complex matrix pair $(F,G)$, according to a given matrix $J$ defining the hyperbolic scalar product. If $J = I$, then these two phases compute the GSVD in parallel very accurately and efficiently.

Figures

Tables

Table 1. Table 1: The datasets under test. For A datasets N L = 121 subscript 𝑁 𝐿 121 N_{L}=121 , N A = 108 subscript 𝑁 𝐴 108 N_{A}=108 , and m = 2 ​ N L ​ N A = 26136 𝑚 2 subscript 𝑁 𝐿 subscript 𝑁 𝐴 26136 m=2N_{L}N_{A}=26136 , while for B datasets N L = 49 subscript 𝑁 𝐿 49 N_{L}=49 , N A = 512 subscript 𝑁 𝐴 512 N_{A}=512 , and m = 2 ​ N L ​ N A = 50176 𝑚 2 subscript 𝑁 𝐿 subscript 𝑁 𝐴 50176 m=2N_{L}N_{A}=50176 . Also, n = N G 𝑛 subscript 𝑁 𝐺 n=N_{G} .

ID	AuAg	$n$
A1	2.5	03275
A2	3.0	05638

Table 2. Table 2: The average per-atom wall execution time (wtime) of Phase 1 with 32 and 64 threads. Since the routine weights are rounded to the nearest per mil, their sum may not yield 100%. The first weight corresponds to ZHEBPJ , the second one to ZGEMM , and the third one to ZDSCAL s step of Algorithm  1 .

ID	average wtime [s] per atom		routine weights %:%:%
	32 threads	64 threads	32 threads	64 threads
A1	$0.243186$	$0.277876$	$71.9:26.2:1.9$	$66.1:32.1:1.8$
A2	$0.279318$	$0.312057$	$59.6:36.6:3.8$	$53.8:41.4:4.9$
A3	$0.345113$	$0.409165$	$47.9:46.4:5.7$	$41.4:46.2:12.4$
A4	$0.436803$	$0.536776$	$37.8:53.5:8.8$	$31.7:50.9:17.5$
B1	$0.023099$	$0.027365$	$56.3:38.4:5.3$	$49.9:45.6:4.5$
B2	$0.030430$	$0.033248$	$40.6:52.3:7.2$	$38.6:54.8:6.7$
B3	$0.045586$	$0.060505$	$25.6:65.1:9.3$	$20.2:71.7:8.1$
B4	$0.070669$	$0.139536$	$16.1:71.6:12.3$	$8.3:80.9:10.8$

Table 3. Table 3: The wall execution time (wtime) in seconds of the Phase 2 steps: the hyperbolic QR (JQR), the tall-and-skinny QR ( ZGEQR from LAPACK), and the prepermutation of the columns of G ~ ~ 𝐺 \widetilde{G} by P 2 subscript 𝑃 2 P_{2} (max. wtime with 32, 64, and 48 threads). In the last column is the number of 2 × 2 2 2 2\times 2 pivots chosen.

ID	JQR wtime [s]			TSQR wtime [s]			$\tilde{G}{P}_2$	$2\times 2$
	32 thr.	64 thr.	48 thr.	32 thr.	64 thr.	48 thr.	wtime [s]	pivots
A1	$162.52$	$179.02$	$66.52$	$10.53$	$7.19$	$1.80$	$0.31$	$11$
A2	$449.13$	$482.63$	$184.46$	$14.06$	$8.26$	$4.07$	$0.87$	$9$
A3	$1041.92$	$1104.48$	$383.60$	$22.01$	$14.71$	$8.03$	$1.14$	$7$
A4	$2137.63$	$2249.65$	$753.71$	$42.82$	$27.38$	$16.54$	$2.37$	$6$
B1	$181.13$	$177.44$	$64.72$	$12.00$	$12.21$	$3.41$	$0.06$	$16$
B2	$520.15$	$495.70$	$186.48$	$15.63$	$10.29$	$4.28$	$0.69$	$16$
B3	$1263.68$	$1182.83$	$420.02$	$36.66$	$20.22$	$10.29$	$1.59$	$20$
B4	$2780.19$	$2439.50$	$841.83$	$49.04$	$33.12$	$19.35$	$3.41$	$16$

Table 4. Table 4: The wall execution time (wtime), with 32, 64, and 48 threads, of the three algorithm candidates.

ID	Level 1 (VP) wtime [s]		Level 2 (FB) wtime [s]		Level 2 (BO) wtime [s]
	32 thr.	64 thr.	32 thr.	64 thr.	32 thr.	64 thr.	48 thr.
A1	$828.58$	$836.50$	$372.00$	$222.81$	$125.38$	$128.51$	$39.39$
A2	$4528.13$	$4533.41$	$2215.07$	$914.32$	$634.72$	$481.92$	$169.97$
A3	$19094.72$	$19357.67$	$9253.03$	$3539.32$	$2392.01$	$1474.61$	$569.85$
A4	$58910.73$	$54938.04$	$30925.99$	$11449.48$	$8003.41$	$4528.69$	$2128.77$
B1	$254.69$	$266.73$	$135.87$	$101.21$	$46.36$	$50.43$	$16.91$
B2	$1379.01$	$1461.47$	$619.19$	$328.24$	$188.76$	$173.86$	$53.28$
B3	$6377.36$	$6249.12$	$2870.85$	$1142.89$	$789.13$	$529.44$	$190.45$
B4	$22002.45$	$22421.95$	$9620.59$	$4015.04$	$2361.69$	$1506.33$	$579.38$

Table 5. Table 5: The wall execution time (wtime) with 32, 64, and 48 threads, the speedups, and the numbers of block sweeps of Phase 3 on the tall-and-skinny and the square inputs (obtained by Phase 2).

ID	BO tall-and-skinny $(\diamond )$		Phase 2 & BO square $(\star )$			speedup $(\diamond )/(\star )$ & sweeps
	32 thr.	64 thr.	32 thr.	64 thr.	48 thr.	32 thr.	64 thr.
A1	$590.92$	$1035.61$	$298.71$	$315.04$	$107.73$	$1.978;15|14$	$3.287;15|14$
A2	$1813.66$	$1600.17$	$1098.79$	$973.41$	$358.56$	$1.651;17|16$	$1.644;16|16$
A3	$5393.96$	$3457.28$	$3457.09$	$2594.89$	$961.61$	$1.560;19|17$	$1.332;17|16$
A4	$13293.75$	$7734.08$	$10185.98$	$6808.09$	$2899.21$	$1.305;19|18$	$1.136;18|18$
B1	$609.11$	$992.56$	$239.55$	$240.14$	$85.10$	$2.543;12|13$	$4.133;12|12$
B2	$1665.46$	$1716.25$	$725.15$	$680.55$	$244.13$	$2.297;14|13$	$2.522;13|13$
B3	$3595.73$	$3126.35$	$2091.06$	$1733.91$	$620.94$	$1.720;15|15$	$1.803;14|14$
B4	$8031.08$	$6147.36$	$5194.34$	$3981.79$	$1440.84$	$1.546;15|15$	$1.544;15|15$

Table 6. Table 6: The wall times for the explicit formation of H 𝐻 H and S 𝑆 S , combined with those for a LAPACK’s generalized Hermitian eigensolver ( ZHEGV or ZHEGVD ), and the speedup vs.  ( ⋆ ) ⋆ (\star) , with 64 and 48 threads.

ID	# of	(max. of 2 runs for $H,S$) wall time [s] for				total wall time [s] with		speedup vs. $(\star )$
	thrs.	$H={\tilde{F}}^{\ast }\tilde{J}\tilde{F}$	$S={\tilde{G}}^{\ast }\tilde{G}$	ZHEGV	ZHEGVD	ZHEGV $(◁)$	ZHEGVD $(▷)$	$(\star )/(◁)$	$(\star )/(▷)$
A1	64	$1.490$	$1.519$	$11.896$	$5.573$	$14.905$	$8.317$	$21.137$	$37.881$
	48	$1.316$	$0.536$	$2.875$	$1.800$	$4.727$	$3.185$	$22.793$	$33.825$
A2	64	$3.959$	$3.323$	$38.197$	$18.890$	$45.479$	$25.913$	$21.403$	$37.565$
	48	$2.368$	$1.277$	$12.722$	$8.504$	$16.349$	$12.148$	$21.931$	$29.515$
A3	64	$11.642$	$7.696$	$100.320$	$56.839$	$119.448$	$75.966$	$21.724$	$34.158$
	48	$5.808$	$3.201$	$53.464$	$37.089$	$62.472$	$46.098$	$15.393$	$20.860$
A4	64	$24.399$	$15.369$	$284.750$	$168.956$	$324.281$	$208.513$	$20.994$	$32.651$
	48	$12.692$	$7.095$	$160.374$	$119.576$	$180.157$	$139.290$	$16.093$	$20.814$
B1	64	$1.348$	$1.909$	$5.399$	$2.573$	$8.656$	$5.575$	$27.742$	$43.076$
	48	$0.866$	$0.507$	$1.055$	$0.603$	$2.427$	$1.954$	$35.058$	$43.550$
B2	64	$4.401$	$3.786$	$17.015$	$8.750$	$24.987$	$16.648$	$27.236$	$40.879$
	48	$2.270$	$1.281$	$4.284$	$2.791$	$7.799$	$6.342$	$31.302$	$38.493$
B3	64	$11.213$	$10.047$	$45.139$	$24.654$	$63.692$	$45.914$	$27.223$	$37.764$
	48	$5.508$	$3.104$	$15.564$	$12.399$	$24.162$	$20.904$	$25.699$	$29.704$
B4	64	$23.092$	$16.446$	$114.787$	$62.469$	$150.798$	$102.006$	$26.405$	$39.035$
	48	$12.003$	$6.790$	$51.817$	$41.265$	$70.570$	$59.836$	$20.417$	$24.080$

Table 7. Table 7: The wall execution time (wtime) and the speedup of the Phases 3 (with J = I 𝐽 𝐼 J=I ) and 4 versus the ZGGSVD3 LAPACK routine on the set C, with 64 threads and n 𝑛 n denoting the order of the matrices.

$n$	ZGGSVD3	Phase 3	Phase 4	Phases 3 & 4	speedup
	wtime [s] $(\bullet )$	wtime [s] & sweeps	wtime [s]	wtime [s] $(\circ )$	$(\bullet )/(\circ )$
$1000$	$399.47$	$8.32;13$	$2.84$	$11.16$	$35.78$
$2000$	$5935.03$	$35.52;14$	$20.28$	$55.80$	$106.37$
$3000$	$21880.11$	$100.71;16$	$69.48$	$170.19$	$128.57$
$4000$	$54233.01$	$191.32;16$	$180.70$	$372.02$	$145.79$
$5000$	$107424.26$	$332.29;17$	$327.97$	$660.26$	$162.70$

Equations

$$H$$

$$S$$

$$H_{AA}=\sum _{a=1}^{N_A}{A}_a^{\ast }{T}_a^{[AA]}{A}_a,H_{AB+BA+BB}=\sum _{a=1}^{N_A}\left(B_a^{\ast }{Z}_a+Z_a^{\ast }{B}_a\right),$$

$${\hat{H}}_{\mathrm{KS}}{\psi }_i(\text{r})=\left[-\frac{{\hslash }^2}{2m_e}{\nabla }_{\text{r}}^2+V[n(r)](\text{r})\right]{\psi }_i(\text{r})={ϵ}_i{\psi }_i(\text{r}),i=1,\dots ,N_e.$$

$$n(r)=\sum _{i=1}^{N_e}|{\psi }_i(\text{r}){|}^2.$$

$${\psi }_i(\text{r})=\sum _{t=1}^{N_G}{c}_{t,i}{\phi }_t(\text{r}).$$

$${\phi }_t(\text{r})=\{\begin{array}{ll}\sum _{l=0}^{l_{\max }}\sum _{m=-l}^l\left[A_{(l,m),a,t}{u}_{l,a}(r)+B_{(l,m),a,t}{\dot{u}}_{l,a}(r)\right]Y_{l,m}({\hat{\mathbf{r}}}_a),& a^{\mathrm{th}}\text{MT}\\ \frac{1}{\sqrt{\Omega }}\exp (\mathrm{i}{\mathbf{k}}_t\cdot \text{r}),& \text{INT}.\end{array}$$

$$\sum _{t=1}^{N_G}(H{)}_{t^{\prime },t}{c}_{t,i}={ϵ}_i\sum _{t=1}^{N_G}(S{)}_{t^{\prime },t}{c}_{t,i}.$$

$$(H{)}_{t^{\prime },t}=\sum _a\iint {\phi }_{t^{\prime }}^{\ast }(\text{r}){\hat{H}}_{\mathrm{KS}}{\phi }_t(\text{r})\mathrm{d}\text{r},(S{)}_{t^{\prime },t}=\sum _a\iint {\phi }_{t^{\prime }}^{\ast }(\text{r}){\phi }_t(\text{r})\mathrm{d}\text{r}.$$

$$(H{)}_{t^{\prime },t}$$

$$+\left(B_{L^{\prime },a,t^{\prime }}^{\ast }{T}_{L^{\prime },L;a}^{[BA]}{A}_{L,a,t}\right)+\left(B_{L^{\prime },a,t^{\prime }}^{\ast }{T}_{L^{\prime },L;a}^{[BB]}{B}_{L,a,t}\right),$$

$$(S{)}_{t^{\prime },t}$$

$$F=U{\Sigma }_{F}X,G=V{\Sigma }_{G}X,{\Sigma }_F\in {\mathbb{R}}^{m\times n},{\Sigma }_G\in {\mathbb{R}}^{p\times n}.$$

$$H$$

$$S$$

$$HZ=SZ\Lambda ;Z:=X^{-1},\Lambda :={\Lambda }_G^{-1}{\Lambda }_F.$$

$$T_a=\left[\begin{array}{cc}{T}_a^{[AA]}& T_a^{[AB]}\\ T_a^{[BA]}& T_a^{[BB]}\end{array}\right]$$

$$H=\sum _{a=1}^{N_A}{H}_a^{\ast }{T}_a{H}_a,S=\sum _{a=1}^{N_A}{S}_a^{\ast }{S}_a,H_a=\left[\begin{array}{c}{A}_a\\ B_a\end{array}\right],S_a=\left[\begin{array}{c}{A}_a\\ U_a{B}_a\end{array}\right].$$

$$H$$

$$S$$

$$T_a=P_a^T{M}_a^{\ast }{D}_a{M}_a{P}_a,$$

$$H={\tilde{F}}^{\ast }\tilde{J}\tilde{F},{\tilde{F}}^{\ast }=[{\tilde{H}}_1^{\ast },\dots ,{\tilde{H}}_{N_A}^{\ast }],\tilde{J}=\mathrm{diag}({\tilde{J}}_1,\dots ,{\tilde{J}}_{N_A}).$$

$$P_1\tilde{F}{P}_2=Q_{F}F;{\tilde{Q}}_F^{\ast }\tilde{J}{\tilde{Q}}_F=J,{\tilde{Q}}_F:=P_1^T{Q}_F,$$

$$P_2^T{\tilde{F}}^{\ast }\tilde{J}\tilde{F}{P}_2=F^{\ast }{Q}_F^{\ast }{P}_1\tilde{J}{P}_1^T{Q}_{F}F=F^{\ast }{\tilde{Q}}_F^{\ast }\tilde{J}{\tilde{Q}}_{F}F=F^{\ast }JF.$$

$$P_3(\tilde{G}{P}_2)=Q_{G}G;{\tilde{Q}}_G^{\ast }{\tilde{Q}}_G=I_n,{\tilde{Q}}_G:=P_3^T{Q}_G,$$

$$P_4(\tilde{G}{P}_2)P_5=Q_{G^{\prime }}{G}^{\prime };{\tilde{Q}}_{G^{\prime }}^{\ast }{\tilde{Q}}_{G^{\prime }}=I_n,{\tilde{Q}}_{G^{\prime }}:=P_4^T{Q}_{G^{\prime }},$$

$$\mathrm{Re}({\Sigma }_j)$$

$$\mathrm{Im}({\Sigma }_j)$$

$${\Sigma }_j^{\prime }={\Sigma }_j^{\prime }+{\tilde{J}}_i(\mathrm{Re}(f_i{)}^2+\mathrm{Im}(f_i{)}^2).$$

$$H(w)=I-2w(w^{\ast }{\tilde{J}}_{k}w{)}^{+}{w}^{\ast }{\tilde{J}}_k,$$

$${\tilde{f}}_1^{\ast }{\tilde{J}}_k{\tilde{f}}_1$$

$${\tilde{f}}_1^{\ast }{\tilde{J}}_k{f}_1$$

$${\overline{\stackrel{~}{f}}}_{11}{\tilde{ȷ}}_{11}{c}_1={\bar{c}}_1{\tilde{ȷ}}_{11}{\tilde{f}}_{11},$$

$$w^{\ast }{\tilde{J}}_{k}w$$

Full text

{AMS}

65F15, 65F25, 65Y05, 65Z05

1 Introduction

Density Functional Theory (DFT) is the Standard Model at the base of simulations in condensed matter physics. At the center of most DFT simulations lays the initialization of the Hamiltonian matrix $H$ and its diagonalization. In many DFT methods the form and size of the Hamiltonian depends on the choice of the set of basis functions used to expand the atomic orbitals. When such a basis set is not orthonormal, a Hermitian positive definite overlap matrix $S$ has to be computed and diagonalized simultaneously with $H$; this pair of matrices $(H,S)$ define a generalized Hermitian eigenproblem (or eigenpencil, in short). In a subset of all DFT methods labeled as LAPW, the entries of both $H$ and $S$ are represented as multiple sums and products of smaller matrices with specific properties. We show how to exploit this peculiar representation to solve the generalized eigenvalue problem without explicitly assembling the $H$ and $S$ matrices. Our alternative method solves the eigenpencil using a cascade of phases ending with the Hari–Zimmermann algorithm for a generalized hyperbolic SVD. We demonstrate the scalability of a shared memory version of this method on a number of test cases extracted from concrete DFT simulations. If the matrix $S$ is ill-conditioned, our method has the additional benefit of providing enhanced accuracy while avoiding the failure-prone Cholesky factorization.

The birth of DFT is marked by two fundamental articles by the Nobel prize winner Walter Kohn and his collaborator Lu J. Sham and Pierre Hohenberg [12, 15]. DFT provides an approach to the theory of electronic structure that is alternative to the solution of the Schrödinger equation. While in the latter the emphasis is on a many-electron wave function describing the dynamics of electrons in a multi-atomic system, in DFT the electron density distribution $\operatorname{\mathit{n}}(\hbox{{r}})$ plays a central role. Besides providing a complementary perspective, DFT has made possible the simulation of much larger systems than the conventional multi-particle wave function methods. Depending on the specific DFT method, computing complexity scales at most with the cube of the number of atoms, with ongoing progress towards bringing it down to linear scaling.

Despite being a general theory, DFT can be realized in as many flavors as are the sets of basis functions one can choose from. Two widely spread classes of basis functions build on the simplicity of plane waves to build more complex and rich sets of basis functions, namely Projected Plane Waves (PAW) [23] and Linearized Augmented Plane Waves (LAPW) [26]. The complexity of these sets lays in that they are made up of non-orthogonal basis functions. In the case of LAPW, the set of functions is also overcomplete. The consequence of non-orthogonality is that the matrix $S$, whose entries are the scalar products among all the basis functions of any given finite size set, is usually dense. In the particular case of LAPW methods such matrix is positive definite but could have few singular values quite close to zero. This potential problem is due to the overcompleteness of the basis set and tends to worsen as the number of atoms increases since the number of basis functions grows linearly with it.

In DFT methods the dynamics of the quantum systems is described by a Hamiltonian operator. In practice, the Hamiltonian is translated into a Hermitian matrix $H$ whose size and structure depends on the specific DFT method. This is because the matrix $H$ is the result of the projection of the Hamiltonian operator over the finite set of basis functions of the given method. In the LAPW method, the mathematical form of the functions leads to an expression for both $H$ and $S$ in terms of a sum of smaller matrices over all possible atoms $N_{A}$,

[TABLE]

where $A_{a},B_{a}\in\mathbb{C}^{N_{L}\times N_{G}}$, with $N_{G}$ and $N_{L}$ ($N_{G}\geq N_{L}$) being the size of the basis set and the total number of angular momentum states, respectively. The remaining matrices in (1) are complex, square, of order $N_{L}$, with some additional properties. Matrices $U_{a}$ are real and diagonal, for all $T_{a}^{[AB]}$ holds $(T_{a}^{[AB]})^{\ast}=T_{a}^{[BA]}$, while $T_{a}^{[AA]}$ and $T_{a}^{[BB]}$ are Hermitian. Except for $U_{a}$, the other matrices are in general dense, and can have a range of sizes dictated by the constants $N_{A}$, $N_{G}$, and $N_{L}$ (see section 2 for some their typical range). Despite the formulation above could lend itself to computation through specialized middleware libraries such as the Basic Linear Algebra Subprograms (BLAS), the standard approach followed by most code developers was one based on minimizing memory footprint and FLOP count [5, 16].

Recently, an alternative method for the assembly of the matrices $H$ and $S$ was presented in [8] and further developed in [6]. In their work [8], Di Napoli et al. consolidate the underlying matrix structure of the operations and proceed to encapsulate them in terms of the level 3 BLAS kernels. For instance, to maximize the arithmetic intensity of the computation, matrix $H$ is written as $H=H_{AA}+H_{AB+BA+BB}$, with

[TABLE]

where $Z_{a}=T_{a}^{[BA]}A_{a}+\nicefrac{{1}}{{2}}T_{a}^{[BB]}B_{a}$.

Each of the $Z_{a}$ and $B_{a}$ matrices is then packed in memory in two consecutive 2-dimensional arrays $Z_{*}$ and $B_{*}$, respectively. In the end, the sum $H_{AB+BA+BB}$ is computed by just two ZHER2K BLAS subroutines. A similar procedure holds for the matrix $S$. Once assembled the algebraic dense generalized eigenproblem is solved by standard methods. A Cholesky factorization $LL^{\ast}=S$ is used to reduce the problem to standard form $A\leftarrow L^{-1}AL^{-\ast}$. In turn, the standard problem is solved by a dense direct algorithm such as MRRR [7] provided by the LAPACK library [1], or an iterative eigensolver specialized for DFT computation (e.g., the ChASE library [29]). When the assembled $S$ matrix is ill-conditioned, as it may happen for quantum systems with a large number of atoms ($>100$), the Cholesky factorization may fail and makes it practically very hard to solve the corresponding generalized eigenproblem. This issue is typically solved by the practitioners by modifying the mathematical model so to avoid increasing the basis set, which is the source of an ill-condioned $S$, but thus compromising on the robustness of the DFT approach.

In this work, we propose an numerical method alternative to the physics-based approach for solving the eigenpencil $(H,S)$ without forming the matrices explicitly. The core of the method is based on the generalized hyperbolic singular value decomposition (GHSVD) [2]. Not only such a method solves for the eigenproblem directly without assembling $H$ and $S$, but also could give more accurate results when $S$ is nearly singular. This is possible since the GHSVD decomposition acts directly on the multiplying factors making up $S$, conceivably reducing the singularity down to the square root of the condition number of $S$. As a surplus, if $J$, the matrix of the hyperbolic scalar product, is equal to the identity, the GHSVD reduces to the generalized SVD (GSVD), which is computed very efficiently in parallel.

The paper is subdivided into eight sections. In section 2, we present in more detail the physics of the problem and the mathematical model leading to the expression (1) for $H$ and $S$. Section 3 is devoted to formulating the problem in precise algebraic terms, and provides an overview of our algorithm. The next four sections deal with the four phases of the algorithm, where the first three of them belong to the algorithm proper, and the fourth one completes the computation of the GHSVD and is unrelated to the underlying mathematical physics origin of the problem. Since each phase is an algorithm and a reusable software contribution in its own right, at the end of each section we present the numerical results and the parallelization techniques applied. The paper concludes with a note on the related future work in section 8.

2 The $H$ and $S$ matrices in LAPW methods

At the core of DFT are a set of equations, called Kohn–Sham equations, that have to be solved for each of the single particle wave function $\psi_{i}$

[TABLE]

The peculiarity of these equations is that the Hamiltonian operator $\hat{H}_{\rm KS}$ depends implicitly on all the $\psi$ through the charge density function $\operatorname{\mathit{n}}(\hbox{{r}})$, which makes the entire set of Kohn–Sham equation strongly coupled and non-linear. In particular the function $\operatorname{\mathit{n}}(\hbox{{r}})$ is the sum of the squares of all $\psi$ up to the total number of electrons $N_{e}$ in any given quantum system

[TABLE]

Because the equations (2) are non-linearly coupled, they can be solved only self-consistently: one starts from a reasonable guess for the charge density $\operatorname{\mathit{n}}(\hbox{{r}})_{\rm start}$, computes the potential $V[\operatorname{\mathit{n}}(\hbox{{r}})]$, and solves (2). The resulting functions $\psi_{i}$ and values $\epsilon_{i}$ are then used to compute a new density as in (3), which is compared to the starting one. If the two densities do not match, the self-consistent loop is repeated with a new mixed charge density. The loop stops only when the new and the old density agree up to some defined constant.

So far we have described the general setup. There are many methods that translate this setup into an algorithm, and this is where the various “flavors” of DFT differ. The first difference is in the choice of the set of functions $\varphi_{t}$ used to expand every one particle wave function $\psi_{i}$

[TABLE]

In the LAPW method [13, 28], the configuration space where the atomic cells are defined is divided in two disjoint areas where the wave functions have distinct symmetries: close to the atomic nuclei, solutions tend to be spherically symmetric and strongly varying, while further away from the nuclei, they can be approximated as uniformly oscillating. The qualitative structure of the solution leads to a space composed of non-overlapping spheres—called muffin tins (MT)—separated by interstitial (INT) areas. The complete set of basis functions $\varphi_{t}$ are given by a piece-wise definition for each of the $N_{A}$ atoms and relative surrounding regions.

[TABLE]

In the MT spheres, each basis function depends on specialized radial functions $u_{l,a}$, their derivatives $\dot{u}_{l,a}$ and the spherical harmonics $Y_{l,m}$; the former only depend on the distance r from the MT center, while the latter form a complete basis on the unit sphere defined by $\hat{\bf r}=\nicefrac{{\textbf{r}}}{{|\textbf{r}|}}$ and so depends solely on the MT spherical angles. Despite being piece-wise functions, $\varphi_{t}$ must be continuous and differentiable for each index $t$ and each atomic index $a$. The coefficients $A_{l,m,a},B_{l,m,a}\in\mathbb{C}$ are set to guarantee that $\varphi_{t}\in C^{1}$ for each of the values of the indices $L\equiv(l,m)$ and $a$. The variable $t$ ranges over the size of the plane wave functions set in INT, and is used to label the vector $\hbox{{k}}_{t}$ living in the space reciprocal to r. As such, the momentum $\hbox{{k}}_{t}$ characterizes the specific wave function entering in the basis set. The total size of the basis set is determined by setting a cutoff value ${\bf K}_{\max}\geq\hbox{{k}}_{t}$.

When one substitutes the expansion of $\psi_{i}$ (4) in (2), the Kohn–Sham equations become an algebraic generalized eigenvalue problem that needs to be solved for the $N_{G}$-tuples of coefficients $c_{i}=(c_{1,i},\ldots,c_{N_{G},i})^{T}$

[TABLE]

The complexity of the LAPW basis set is transferred to the definition of the entries of the Hamiltonian and overlap matrices, respectively $H$ and $S$, given by

[TABLE]

By substituting explicitly the functions $\varphi_{t}$ of equation (5) and computing the integrals, one ends up with the following expressions for $H$ and $S$:

[TABLE]

The new matrices $T_{L^{\prime}\!,\mkern 1.0muL;a}^{[\dots]}\in\mathbb{C}^{N_{L}\times N_{L}}$ are dense and their computation involves multiple integrals between the radial basis functions $u_{l,a}$ and the non-spherical part of the potential $V$ multiplied by Gaunt coefficients (for details see [16][26, Ch. 5][8, App.]). As can be seen by simple inspection, equations (6)–(7) are equivalent to equations (1): while the former are written with all indices explicit, the latter have a subset of them implicit which highlights their matrix form.

We conclude with a small excursus on the structure of the self-consistent loop and its computational cost. In the first step, a starting charge density $\operatorname{\mathit{n}}(\hbox{{r}})_{\rm start}$ is used to compute the Kohn–Sham Hamiltonian $H_{\rm KS}$. In a second step the set of basis functions is set up and the set of $A,B$ coefficients is derived. Then, the Hamiltonian $H$ and overlap $S$ matrices are initialized, followed by the fourth step when the generalized eigenvalue problems $Hc=\epsilon\ Sc$ is solved numerically to return the eigenpairs $(C,{\rm diag}(\epsilon))$. Finally a new charge density $\operatorname{\mathit{n}}(\hbox{{r}})$ is computed and convergence is checked before starting a new loop. Out of all the steps above, initializing $H$ and $S$ and solving the eigenproblem accounts for more than 80% of CPU time. Having cubic complexity $\mathcal{O}(N_{G}^{3})$, the eigenproblem solution is usually considered the most expensive of the two. It turns out that generating the matrices may be as expensive. If $N_{A}$ and $N_{L}$, respectively, are the range of the summations $\sum_{a}$ and $\sum_{L}$, then it can be shown that equations (7) and (6) have complexity $\mathcal{O}(N_{A}\cdot N_{L}\cdot N_{G}^{2})$ and $\mathcal{O}(N_{A}\cdot N_{L}\cdot N_{G}\cdot(N_{L}+N_{G}))$. A typical simulation uses approximately $N_{G}$ basis functions, with $N_{G}$ ranging from about $50\cdot N_{A}$ to about $80\cdot N_{A}$, and an angular momentum $l_{\mathrm{max}}\leq 10$, which results in $N_{L}=(l_{\max}+1)^{2}\leq 121$. It follows that the factor $N_{A}\cdot N_{L}$ is roughly of the same order of magnitude as $N_{G}$ so that the generation of $H$ and $S$ also displays cubic complexity $\mathcal{O}(N_{G}^{3})$. In practice, the constants above have values in the following orders of magnitude: $N_{A}=\mathcal{O}(100)$, $N_{G}=\mathcal{O}(1000)$–$\mathcal{O}(10000)$, and $N_{L}=\mathcal{O}(100)$.

3 Problem formulation

Our intention is to keep the matrices $H$ and $S$ in their factored form given in (1). The core of the process, Phase 3, is a one-sided Jacobi-like method for the implicit diagonalization, that computes a hyperbolic analog of the generalized SVD.

Definition 3.1**.**

For the given matrices $F\in\mathbb{C}^{m\times n}$, $m\geq n$, $J\in\mathbb{R}^{m\times m}$, $J=\operatorname{diag}(\pm 1)$, and $G\in\mathbb{C}^{p\times n}$, where $G$ is of full column rank, there exist a $J$-unitary matrix $U\in\mathbb{C}^{m\times m}$ (i.e., $U^{\ast}JU=J$), a unitary matrix $V\in\mathbb{C}^{p\times p}$, and a nonsingular matrix $X\in\mathbb{C}^{n\times n}$, such that

[TABLE]

*The elements of $\Sigma_{F}$ and $\Sigma_{G}$ are zeros, except for the diagonal entries, which are real and nonnegative. Furthermore, $\Sigma_{F}$ and $\Sigma_{G}$ satisfy $\Sigma_{F}^{T}\Sigma_{F}+\Sigma_{G}^{T}\Sigma_{G}=I$. The ratios $\Sigma_{ii}\mathrel{\mathop{:}}=(\Sigma_{F})_{ii}/(\Sigma_{G})_{ii}$ are called the generalized hyperbolic singular values of the pair $(F,G)$. If the pair $(F,G)$ is real, then all matrices in (8) are real. *

We choose to define the generalized hyperbolic SVD (GHSVD) only if the matrix $G$ is of full column rank. This implies $p\geq n$, and there is no need to mention this in the definition. In the case of full column rank $G$, the matrix $S\mathrel{\mathop{:}}=G^{\ast}G$ is positive definite, and the matrix pair $(H,S)$, where $H\mathrel{\mathop{:}}=F^{\ast}JF$, is Hermitian and definite, so it can be simultaneously diagonalized by congruences (see, for example, [22]).

If the GHSVD is computed as in (8), then the generalized eigenvalues and eigenvectors of $(H,S)$ are easily retrieved, since

[TABLE]

Substituting $X^{\ast}=SX^{-1}\Lambda_{G}^{-1}$ in the expression for $H$ above, we get

[TABLE]

Thus, from (9), the generalized eigenvalues $\operatorname{diag}(\Lambda)$ of the matrix pair are the squared generalized hyperbolic singular values, with the signs taken from the corresponding diagonal elements in $J$, i.e., $\operatorname{diag}(\Sigma^{T}J\Sigma)$, and the matrix of the generalized eigenvectors $Z$ is the inverse of the matrix $X$ of the right generalized singular vectors. For theoretical purposes it can be assumed that $\operatorname{diag}(\Lambda)$ is sorted descendingly, though for simplicity it is not the case in our implementation.

An approach that uses the SVD on a matrix factor, instead of the eigendecomposition on the multiplied factors, usually computes small eigenvalues more accurately.

In the first phase of the algorithm, we transform the initial problem by assembling the Hermitian matrices $T_{a}$,

[TABLE]

and factoring them into a form suitable for the GHSVD computation.

After that, we are left with two tall matrices, which have, in our test examples, between $2$ and $22$ times more rows than columns. Since the Jacobi-like SVD algorithms are more efficient if the factors are square, in the second (optional) phase we could preprocess the factors: $F$ by the hyperbolic QR factorization (see [24]), and $G$ by the tall-and-skinny QR factorization (ZGEQR routine from LAPACK), into square ones.

The third phase is a complex version of the implicit Hari–Zimmermann method—a modification of the one-sided real method presented in [21]. The complex transformations, for the two-sided method, were derived by Vjeran Hari in his PhD thesis [9].

3.1 Overview of the algorithm

The sequence of phases of our algorithm is:

The problem is expressed as $H=F_{0}^{\ast}\operatorname{diag}(T_{1},\ldots,T_{N_{A}})F_{0}$, $S=\widetilde{G}^{\ast}\widetilde{G}$, the matrices $F_{0}$ and $\widetilde{G}$ are assembled, and the matrices $T_{a}$, formed from $T_{a}^{[AA]}$, $T_{a}^{[AB]}$, $T_{a}^{[BA]}$, and $T_{a}^{[BB]}$, are simultaneously factored by the Hermitian indefinite factorization with complete pivoting, reformulating $H$ as $H=\widetilde{F}^{\ast}\widetilde{J}\widetilde{F}$, with $\widetilde{J}=\operatorname{diag}(\pm 1)$. 2. 2.

Optionally, the tall-and-skinny matrices $\widetilde{F}$ and $\widetilde{G}$ are shortened: $\widetilde{F}$ by the indefinite, $\widetilde{J}$-QR factorization to obtain the square factor $F$ and a new signature matrix $J$, and $\widetilde{G}$ by the QR factorization to obtain the square factor $G$. 3. 3.

The GEVD of $(H,S)$ is computed by the $J$-GHSVD of $(F,G)$ (or the $\widetilde{J}$-GHSVD of $(\widetilde{F},\widetilde{G})$ if the Phase 2 is skipped) by the implicit Hari–Zimmermann method. 4. 4.

Optionally, the GHSVD process is formally completed by explicitly computing the right generalized singular vectors $X$ from the generalized eigenvector matrix $Z$.

3.2 Testing environment and data

The testing environment consists of a node with an Intel Xeon Phi 7210 CPU, running at 1.3 GHz with Turbo Boost turned off, in Quadrant cluster mode with 96 GiB of RAM and 16 GiB of flat-mode MCDRAM, under 64-bit CentOS Linux 7 with the Intel compilers (Fortran, C) and Math Kernel Library (MKL) version 19.0.5.281, and GNU Fortran 8.3.1 for the error testing.

The software code, freely available in https://github.com/venovako/FLAPWxHZ repository, of all the phases presented in this paper is written mostly in Fortran, with some auxiliary parts in C, while the parallelization relies on the OpenMP constructs.

The phases are meant to be run in a sequence, where each phase is executed as a separate process with several OpenMP threads. Since the modern compute nodes generally have enough memory to hold all required data, the algorithms are implemented for the shared memory, but at least the algorithms for the Phases 1, 3, and 4 can be transformed into distributed-memory ones, should the volume of data so require.

In testing it was established that each thread should be bound to its own physical CPU core, with OMP_PROC_BIND=SPREAD placement policy. The double precision and the double-complex BLAS and LAPACK routines were provided by the thread-parallel MKL, but with only one (i.e., the calling) thread allowed per call, except for the ZSWAP, ZROT, and ZGEQR routines in Phase 2, ZGETC2 routine in Phase 4, and ZGEMM, ZHERK, ZHEGV, and ZHEGVD routines in subsection 6.5.4, where the MKL was allowed to use the test’s upper limit on the number of threads. The nested parallelism is therefore possible but not required in our code. Hyper Threading was enabled but not explicitly utilized, though nothing precludes a possibility that on a different architecture the BLAS or LAPACK calls could benefit from some form of intra-core symmetric multi-threading. A refined thread placement policy OMP_PROC_BIND=SPREAD,CLOSE might then allow better reuse of data in the cache levels shared among the threads of a core.

Apart from the maximal number of threads set to the number of CPU cores in a node, the tests were also performed with half that number, to assess the effects on the computational time of the larger block sizes and the availability of the whole L2 data cache (1 MiB, shared among two cores) to a thread. The algorithms do not constrain the number of threads in principle, but are not intended to be used single-threadedly.

Our main test node has 64 cores, but a subset of the tests were repeated on a faster JUWELS [14] node, with two Intel Xeon Platinum 8168 CPUs, running at 2.7 GHz with 1 MiB L2 cache per each of $2\times 24=48$ cores, with a similar software setup, for a comparison of the GHSVD and the generalized eigendecomposition approaches (see subsection 6.5.4). When the results obtained on JUWELS are shown, the test’s number of threads is emphasized (e.g., 48), to distinguish them from the main results.

Another hardware feature targeted is the SIMD vectorization: each core of both machines has a private L1 data cache of 32 kiB with a line size of 64 B, and equally wide (e.g., 8 double precision floating-point numbers) vector registers upon which a subset of AVX-512 instructions is capable to operate in the SIMD fashion. The vectorization is employed both implicitly, by aligning the data to the cache line size whenever possible and instructing the compiler to vectorize the loops, and semi-explicitly, as will be described in the following sections. The code is parametrized by the maximal SIMD length (i.e., the number of 8 B lanes in the widest vector register type) $\mathsf{v}$, and it vectorizes successfully on other architectures (e.g., on AVX2, with $\mathsf{v}=4$).

Under an assumption that the compiler-generated floating-point reductions (e.g., those of the SUM Fortran intrinsic) obey the same order of operations in each run, and due to the alignment enforced as above, the algorithms should be considered conditionally reproducible, in a sense that the multiple runs of the same executables on the same data in the same environment should produce bitwise-identical results.

3.2.1 Datasets

Each dataset under test contained all matrix inputs ($A_{a}$, $B_{a}$, $U_{a}$, $T_{a}^{[AA]}$, $T_{a}^{[BA]}$, $T_{a}^{[BB]}$) for a single problem instance. With 8 datasets from Table 1 we believe to have a representative coverage of the small-to-medium size problems from practice.

As already mentioned in section 2, the maximum value of the momentum ${\bf K}_{\max}$ which appears as an index to the dataset label (e.g., AuAg_2.5) determines the size of the basis functions set $N_{G}$. This is why datasets with same label (e.g., AuAg) but different index (e.g., 2.5 vs. 3.0) have differing values for $N_{G}$. In the following, the datasets are referred to by their IDs.

4 Phase 1 – simultaneous factorizations of $T_{a}$ matrices

The goal of this section is to rewrite the problem (1) in a form suitable for GHSVD computation.

4.1 Problem reformulation

The first step is to write (1) as

[TABLE]

Furthermore, the matrices in (11) can be expressed as

[TABLE]

In (12), $\operatorname{diag}(T_{1},\ldots,T_{N_{A}})$ stands for a block-diagonal matrix with the prescribed diagonal blocks $T_{a}$, $a=1,\ldots,N_{A}$ from (10). Newly defined matrices have the following dimensions: $H_{a},S_{a}\in\mathbb{C}^{(2N_{L})\times N_{G}}$, $T_{a}\in\mathbb{C}^{(2N_{L})\times(2N_{L})}$, $F_{0},\widetilde{G}\in\mathbb{C}^{(2N_{A}N_{L})\times N_{G}}$, and $T\in\mathbb{C}^{(2N_{A}N_{L})\times(2N_{A}N_{L})}$. From now on, let $m\mathrel{\mathop{:}}=2N_{A}N_{L}$, and $n\mathrel{\mathop{:}}=N_{G}$.

To efficiently exploit the structure of the problem, the matrix $T$ needs to be diagonal, with its diagonal elements equal to either $1$ or $-1$ (possibly with some zeros in the case of a singular $T$). There is no theoretical obstacle to apply the simultaneous ($J$-)orthogonalization in the computation of the GHSVD on the matrices $F_{0}$, $\widetilde{G}$, and $T$ implicitly, but the repeated multiplication (in each reduction step) by $T$ is slow. Therefore, $T$ should be either factored concurrently, by using a modified version of the Hermitian indefinite factorization of all $T_{a}$ blocks, or diagonalized concurrently: all factorizations (or diagonalizations) are independent of each other and can proceed in parallel. Since the diagonalization, compared to the Hermitian indefinite factorization, is a slower process, our choice is to factor all the diagonal blocks $T_{a}$.

4.2 Hermitian indefinite factorization

Each $T_{a}$ is factored by the algorithm described in [27]. The algorithm for each $T_{a}$ consists of the Hermitian indefinite factorization with a suitable pivoting [4], followed by the transformation of the block-diagonal matrix. Such factorization has the following form

[TABLE]

where $P_{a}$ is a permutation (in the LAPACK sense), $M_{a}$ is upper triangular, and $D_{a}$ is block-diagonal, with diagonal blocks of order $1$ or $2$.

Then, $D_{a}$ is transformed into $\widehat{J}_{a}=\operatorname{diag}(\pm 1)$. If $D_{a}$ has a diagonal block of order $1$ at position $k$, then $\widehat{J}_{a}$ stores the sign of this block in its $k$th diagonal element, and the $k$th row of $M_{a}$ is scaled by $|(D_{a})_{kk}|^{1/2}$. In the case of a (Hermitian) pivot block of order $2$, this block is diagonalized by a Jacobi rotation $R_{k}$, and the corresponding two rows of $M_{a}$ in (13) are multiplied by $R_{k}$. Two transformations of the new diagonal elements of $D_{a}$ are then performed, as above. To speed-up the process, the rotation and the scaling of two rows of $M_{a}$ are combined and then applied as a single transformation.

The outer permutations $P_{a}$ are generated starting from the identity, and stored as the partial permutations of the principal submatrices, as in LAPACK, according to the pivoting of choice. Since the matrices $T_{a}$ are of a relatively small order, our choice is the complete pivoting from [4].

4.3 Postprocessing

After the factorization, with a postprocessing step we obtain $T_{a}=\widehat{M}_{a}^{\ast}\widehat{J}_{a}\widehat{M}_{a}$, where $\widehat{M}_{a}=M_{a}P_{a}$, and $\widehat{M}_{a}$ does not need to remain triangular.

Finally, by applying an inner permutation $\widehat{P}_{a}$, $\widehat{J}_{a}$ can be rearranged into a diagonal matrix $\widetilde{J}_{a}$, where the positive signs precede the negative ones on the diagonal. This property of $\widetilde{J}_{a}$ matrices can be exploited to speed-up computation of the hyperbolic scalar products $x^{\ast}Jx$ in the subsequent phases (see subsections 4.5.1, 5.1, and 6.4.2). Let the whole factorization routine described thus far be called ZHEBPJ. Then, $T_{a}=\widetilde{M}_{a}^{\ast}\widetilde{J}_{a}\widetilde{M}_{a}$, $\widetilde{M}_{a}=\widehat{P}_{a}\widehat{M}_{a}$, and $H_{a}$ is multiplied by $\widetilde{M}_{a}$ from the left as $\widetilde{H}_{a}=\widetilde{M}_{a}H_{a}$.

After such preprocessing, $H$ from (12) is written as

[TABLE]

In datasets A, each $\widetilde{J}_{a}$ has 3 positive and 239 negative signs. In datasets B, a non-consecutive half of $\widetilde{J_{a}}$ matrices are positive definite, and the others are negative definite.

4.4 Implementation and testing

The computational tasks for different indices $a$ are fully independent, and are performed in parallel such that each thread is responsible for one or more indices $a$, as indicated in the pseudocode of Algorithm 1.

Each thread, in turn, performs the three steps for its index $a$ sequentially, up to a possible usage of a parallel BLAS in the first two steps. The last step, i.e., computing $U_{a}B_{a}$ to assemble $\widetilde{G}$, is a loop with the independent iterations, and could be done in parallel, using the nested parallelism within each thread, should $N_{G}$ be large enough and should also the newly spawned threads for that loop have enough computational resources available to warrant the overhead of the additional thread management.

Each thread is responsible for allocating (MCDRAM is not explicitly used) and accessing the memory for the data it processes, so the data locality is achievable whenever each NUMA node has enough storage. The algorithm is thus viable in the heavily non-uniform memory access settings, like the Intel Xeon Phi’s SNC-4 mode.

In a distributed memory setting (e.g., using the MPI processes), the assembling of $\widetilde{F}$, $\widetilde{J}$, and $\widetilde{G}$ can be done by assigning to each process a (not necessarily contiguous) subrange of the iteration range of the for-all loop from Algorithm 1, while inside the process all atoms assigned to it are processed exactly as above, within an OpenMP parallel-do loop. The matrices $\widetilde{F}$, $\widetilde{J}$, and $\widetilde{G}$ would then end up being distributed in the chunks corresponding to the chosen subranges among the processes.

4.4.1 Testing

In Table 2 the average per-atom wall execution time of Phase 1 is shown.

The results suggest that it is beneficial to have more L2 data cache available per thread, as is the case with 32 threads overall. In the breakdown of the weights (i.e., percentages of time taken) of each computational step it is confirmed that ZGEMM starts to dominate the other computational steps of Algorithm 1 as the ratio $n/m$ increases. It is a strong indication that even a procedure more expensive than ZHEBPJ, such as a diagonalization of $T_{a}$, may be applied on the datasets having a square-like shape, without considerably degrading the relative performance of Phase 1.

For a fully vectorized, cache-friendly alternative to applying ZDSCAL with a non-unit stride in Algorithm 1 please refer to section S.1 of the supplementary material.

4.5 An alternative way forward

After this phase has completed, one can proceed as described in the rest of the paper, should the condition numbers of (the yet unformed) matrices $H$ and $S$ be large enough to severely affect the accuracy of a direct solution of the generalized Hermitian eigenproblem with the pair $(H,S)$.

An alternative and more time-efficient way to proceed would be to explicitly form $H$ and $S$. For $S=\widetilde{G}^{\ast}\widetilde{G}$, one ZHERK call would suffice. For $H=\widetilde{F}^{\ast}\widetilde{J}\widetilde{F}$, a copy of $\widetilde{F}$ should be made, and that copy’s rows should be scaled in parallel by the diagonal elements of $\widetilde{J}$. One ZGEMM call on $\widetilde{F}^{\ast}$ and $\widetilde{J}\widetilde{F}$ then completes the formation of $H$. After that, an efficient solver for the generalized Hermitian eigenproblem can be employed on $(H,S)$, such as ZHEGV or ZHEGVD from LAPACK, as shown in subsetion 6.5.4.

4.5.1 Row scaling

A cache-friendly implementation of the row scaling by $\widetilde{J}$ is to iterate sequentially over the rows of a fixed column $j$, and change the sign of each element $\widetilde{F}_{ij}$ for which $\widetilde{J}_{ii}=-1$, while the outer parallel-do loop iterates over all column indices $j$. However, that implementation can be optimized further.

If $\widetilde{J}$ has its diagonal partitioned into (regularly or irregularly sized) blocks of the same sign, then it can be compactly encoded as a sequence of pairs $(i_{-},l)_{k}$, one for each block of negative signs, where $i_{-}$ is the first index belonging to a block $k$, and $l\geq 1$ is the block’s length. The iteration over all rows and the conditional sign changes as above can be replaced by iteration over all such blocks. For each block, iterate sequentially in the range of indices $i$ from $i_{-}$ to $i_{-}+l-1$, and change the signs unconditionally, thus eliminating the conditional branching based on the sign of $\widetilde{J}_{ii}$.

Such run-length-like encoding is employed in Phase 3, where it also accelerates the hyperbolic dot products in the case where the positive signs precede the negative ones on the diagonal of a sign matrix (i.e., at most one negative block exists) given by ZHEBPJ when forming the square factors for the inner Hari–Zimmermann method.

5 Phase 2 – optional $(J,I)$ URV factorization

The one-sided Jacobi-type algorithms are fastest if they work on square matrices, since the column dot-products and updates are the shortest possible. If the square factors $F$, $G$, and the corresponding $J$ of the matrix pair $(\widetilde{F}^{\ast}\widetilde{J}\widetilde{F},\widetilde{G}^{\ast}\widetilde{G})$, can be found, instead of the rectangular factors $\widetilde{F}$, $\widetilde{G}$, and the corresponding $\widetilde{J}$, we expect that the overhead of such a shortening will be less than the computational time saved by avoiding the rectangular factors. To this end, matrix $\widetilde{F}$ is shortened by using the hyperbolic QR factorization (also called the JQR factorization), according to the given $\widetilde{J}$:

[TABLE]

where $F\in\mathbb{C}^{n\times n}$ is block upper triangular with diagonal blocks of order $1$ or $2$, $J=\operatorname{diag}(\pm 1)\in\mathbb{R}^{n\times n}$ is the shortened signature matrix, $P_{2}\in\mathbb{R}^{n\times n}$ and $P_{1}\in\mathbb{R}^{m\times m}$ are the column and the row permutation matrix, respectively. Our application does not use $Q_{F}\in\mathbb{C}^{m\times n}$, so it is not explicitly formed. From (14) it holds

[TABLE]

Since the JQR requires both row and column pivoting (see [24]), matrix $\widetilde{G}$, with its columns prepermuted according to $P_{2}$, the column pivoting of the JQR, will then be factored by the ordinary (tall-and-skinny) QR factorization (e.g., by the LAPACK routine ZGEQR). The latter QR factorization does not employ column pivoting, but in principle the row pivoting or presorting may be used:

[TABLE]

where $G\in\mathbb{C}^{n\times n}$ is upper triangular, and $Q_{G}\in\mathbb{C}^{m\times n}$, which is not needed in our application. We also do not depend on the special forms of $F$ and $G$ later on.

From (8) it follows that the $\widetilde{J}$-GHSVD of $\widetilde{F}$ and $\widetilde{G}$, and that of $\widetilde{F}P_{2}$ and $\widetilde{G}P_{2}$, differ only in the column permutation of the right singular vectors, i.e., $\widetilde{X}=XP_{2}^{T}$, or, from (9), the row permutation of the eigenvectors, i.e., $\widetilde{Z}=P_{2}Z$, while $\Sigma$, and thus $\Lambda$, stay the same. Therefore, the square factors $F$, $G$, and $J$ can be used in place of $\widetilde{F}$, $\widetilde{G}$, and $\widetilde{J}$ throughout the rest of the computation, and then the results could be easily converted back to the ones of the original problem.

Whenever $\widetilde{G}$ might be badly conditioned, Phase 2 could be skipped, or we should resort to a slower but more stable QR factorization of $\widetilde{G}P_{2}$ with the column pivoting,

[TABLE]

where $P_{5}$ has to be applied back to $F$, and $P_{4}$ comes from an optional row pivoting. The column-pivoted QR factorization is provided by the LAPACK routine ZGEQP3.

Then, $F^{\prime}\mathrel{\mathop{:}}=FP_{5}$ and $G^{\prime}$ (provided that it is not rank-deficient according to a user-defined tolerance) could be substituted for $F$ and $G$ in the rest of the computation. For $X$ (or $Z$) thus obtained it holds $\widetilde{X}P^{\prime}=X$ (or, $P^{\prime}Z=\widetilde{Z}$), where $P^{\prime}\mathrel{\mathop{:}}=P_{2}P_{5}$. Such an approach is not required for our datasets, and therefore it was not tested.

5.1 $\widetilde{J}$-dot products and norms

Since $\widetilde{J}$ can in principle contain the positive and the negative signs in any order, and vectorization is strongly desired, a $\widetilde{J}$-dot product of two vectors, $f^{\ast}\widetilde{J}g$, is computed as $\mathsf{v}$ piecewise sums $\Sigma_{j}$, $1\leq j\leq\mathsf{v}$,

[TABLE]

where $i$ starts with a value of $j$ and increments in steps of $\mathsf{v}$ up to $m$. The components $\operatorname{Re}(\Sigma)$ and $\operatorname{Im}(\Sigma)$ of the resulting $\Sigma$ are obtained by SUM-reducing $\operatorname{Re}(\Sigma_{j})$ and $\operatorname{Im}(\Sigma_{j})$, respectively. Similarly, the square of the $\widetilde{J}$-norm, $f^{\ast}\widetilde{J}f$ is computed by SUM-reducing $\Sigma_{j}^{\prime}$, where

[TABLE]

The square of the $\widetilde{J}$-“norm” of a vector thus obtained can be positive or negative, with a possibility of cancellations inadvertently occurring in the summations. It is an open question how to compute the squares of the $\widetilde{J}$-“norms” both efficiently and accurately, though one possible speed improvement might be to encode $\widetilde{J}$ as described in subsection 4.5.1 and simplify the above three piecewise summations accordingly.

5.2 Pivoting

To achieve the maximal numerical stability, the JQR factorization is usually performed with complete pivoting. In the first step the pivot column(s) are chosen from the $\widetilde{J}$–Grammian matrix $H=\widetilde{F}^{\ast}\widetilde{J}\widetilde{F}$, and later on, in the $k$th step, from the $\widetilde{J}_{k}$–Grammian matrix $H_{k}=\widetilde{F}_{k}^{\ast}\widetilde{J}_{k}\widetilde{F}_{k}$, where $\widetilde{F}_{k}$ the is a part of the matrix yet to be reduced, and $\widetilde{J}_{k}$ is the matrix of signs that corresponds to the unreduced matrix $\widetilde{F}_{k}$ (see Fig. 1). The complete pivoting in the first step needs formation of the whole $H$, i.e., $\mathcal{O}(mn^{2})$ floating-point operations. Such an approach, consistently implemented throughout the algorithm, leads to $\mathcal{O}(m^{2}n^{2})$ operations solely for the choice of pivots. Therefore, we relaxed the pivoting strategy to the diagonal pivoting supplemented with the partial pivoting [3, Algorithm C].

5.2.1 Diagonal and partial pivoting

First, $n-k+1$ squares of the $\widetilde{J}_{k}$-norms $h_{ii}^{[k]}\mathrel{\mathop{:}}=\tilde{f}_{i}^{\ast}\widetilde{J}_{k}\tilde{f}_{i}$, where $\tilde{f}_{i}$ is the $i$th column of $\widetilde{F}_{k}$ and $i\geq 1$, are computed in a parallel-do loop over $i$ and stored in a work array. Since all columns have the length of $m-k+1$, each parallel loop iteration executes (sequentially) in approximately the same time and the work is therefore well balanced among the threads.

Let $j\geq 1$ be the smallest index such that $|h_{jj}^{[k]}|\geq|h_{ii}^{[k]}|$ for all $i$. If $j>1$, the $k$th and the $(j+k-1)$th columns of $\widetilde{F}$ (and thus also the first and the $j$th columns of $\widetilde{F}_{k}$) are swapped. If $k=n$, the column pivoting is completed.

Otherwise, $n-k$ $\widetilde{J}_{k}$-dot products $h_{1j}^{[k]}\mathrel{\mathop{:}}=\tilde{f}_{1}^{\ast}\widetilde{J}_{k}\tilde{f}_{j}$, where $\tilde{f}_{j}$ is the $j$th column of $\widetilde{F}_{k}$ and $j>1$, are computed in a parallel-do loop over $j$ and stored in a complex workspace, while their magnitudes $|h_{1j}^{[k]}|$ are placed in a real workspace. Same as above, this work is well balanced among the threads.

Let $i>1$ be the smallest index such that $|h_{1i}^{[k]}|\geq|h_{1j}^{[k]}|$ for all $j>1$. As in [3], if $|h_{11}^{[k]}|\geq\alpha|h_{1i}^{[k]}|$, with $\alpha\mathrel{\mathop{:}}=(1+\sqrt{17})/8$, the column pivoting in the step $k$ is completed.

Otherwise, $n-k$ $\widetilde{J}_{k}$-dot products $h_{il}^{[k]}\mathrel{\mathop{:}}=\tilde{f}_{i}^{\ast}\widetilde{J}_{k}\tilde{f}_{l}$, where $\tilde{f}_{l}$ is the $l$th column of $\widetilde{F}_{k}$ and $i\neq l\geq 1$, are computed in a parallel-do loop over $l$ and their magnitudes $|h_{il}^{[k]}|$ are stored in a real workspace. This work is only slightly imbalanced among threads, since for $l=i$ a thread assigned to the iteration sets $|h_{ll}^{[k]}|=0$, excluding the value and its index from the search for a maximum, unless all other values are also 0.

Let $j\geq 1$ be the smallest index such that $|h_{ij}^{[k]}|\geq|h_{il}^{[k]}|$ for all $l$. As in [3], if $|h_{11}^{[k]}||h_{ij}^{[k]}|\geq\alpha|h_{1i}^{[k]}|^{2}$, the column pivoting for the step $k$ is completed; else, if $|h_{ii}^{[k]}|\geq\alpha|h_{ij}^{[k]}|$, the $k$th and the $(i+k-1)$th columns of $\widetilde{F}$ (and thus also the first and the $i$th columns of $\widetilde{F}_{k}$) are swapped and the column pivoting in the step $k$ is completed.

Otherwise, a $2\times 2$ pivot is chosen, by taking the first column of $\widetilde{F}_{k}$ and swapping the $(k+1)$th and $(i+k-1)$th columns of $\widetilde{F}$ (and thus also the second and the $i$th columns of $\widetilde{F}_{k}$), if $i\neq 2$; else, the second pivot column is already in place.

The pivot column(s) have thus been brought to the front of the matrix $\widetilde{F}_{k}$ by at most two column swaps. The ensuing row pivoting is explained further below.

5.3 Hyperbolic Householder reflectors

If a single pivot is chosen, the first column $\tilde{f}_{1}$ of $\widetilde{F}_{k}$ is reduced by a hyperbolic Householder reflector [25, Theorem 4.4] to a vector $f_{1}=c_{1}e_{1}$, where $c_{1}\in\mathbb{C}$ and $e_{1}$ is the first vector of the canonical base. A variant of [25, Theorem 4.4] for the hyperbolic scalar product and a simple shape of $f_{1}$ follows.

Theorem 5.1**.**

Let $\widetilde{J}_{k}$ be a hyperbolic scalar product matrix of order $\ell$. Let $\tilde{f}_{1},f_{1}\in\mathbb{C}^{\ell}$ be two distinct vectors. There exists a basic $\widetilde{J}_{k}$ reflector $H(w)$,

[TABLE]

such that $H(w)\tilde{f}_{1}=f_{1}$ if and only if $\tilde{f}_{1}$ and $f_{1}$ satisfy the $\widetilde{J}_{k}$-isometry and $\widetilde{J}_{k}$-symmetry property, respectively

[TABLE]

*and $d=\tilde{f}_{1}-f_{1}\neq 0$ is nondegenerate, i.e., $d^{\ast}\widetilde{J}_{k}d\neq 0$. Furthermore, whenever $H(w)$ exists, it is unique. $H(w)$ can be generated by any $w\in\mathbb{C}^{\ell}$ such that $w=\lambda d$, $\lambda\in\mathbb{C}\setminus\{0\}$. Finally, the same remains valid if we replace $f_{1}$ by $-f_{1}$, and $d$ by $s=f_{1}+\tilde{f}_{1}$.

If $\tilde{f_{1}}=f_{1}$, there is nothing to do in this step, so we take $H(w)=I$. Otherwise, since we want to obtain $f_{1}$ in the form $f_{1}=c_{1}e_{1}$, equation (16) is equivalent to a requirement that the sign of $\tilde{\jmath}_{11}$, the first diagonal element of $\widetilde{J}_{k}$, is equal to the sign of $\tilde{f}_{1}^{\ast}\widetilde{J}_{k}\tilde{f}_{1}=|c_{1}|^{2}\tilde{\jmath}_{11}=f_{1}^{\ast}\widetilde{J}_{k}f_{1}$. As $\tilde{f}_{1}^{\ast}\widetilde{J}_{k}\tilde{f}_{1}$ has already been computed by the diagonal pivoting (see subsection 5.2.1), it is trivial to check if the requirement holds.

If not, it can be shown that there exists at least one element with the correct sign in $\widetilde{J}_{k}$. When there is more than one such element, our implementation sequentially finds the one that corresponds to the largest element in $\tilde{f}_{1}$ by magnitude (say, $l$th). By permuting the diagonal of $\widetilde{J}_{k}$, this element can be brought to the first diagonal position. This implies a corresponding row permutation of $\widetilde{F}_{k}$ that swaps the first and the $l$th rows of $\widetilde{F}_{k}$. For that, we employ the parallel ZSWAP routine, but should the rows be short enough, a sequential version of the routine could be considered instead.

Relation (16) implies that $|c_{1}|=|\tilde{f}_{1}^{\ast}\widetilde{J}_{k}\tilde{f}_{1}|^{1/2}$. In general, $c_{1}$ is a complex number, $c_{1}=re^{\mathrm{i}\delta}$, and we have already determined $r=|c_{1}|$. It remains to find $\delta=\arg(c_{1})$.

From (17) it follows

[TABLE]

where $\tilde{f}_{11}\mathrel{\mathop{:}}=r_{1}e^{\mathrm{i}\delta_{1}}$ is the first element in $\tilde{f}_{1}$. Relation (18) can be divided by $\tilde{\jmath}_{11}$ and written as $r_{1}re^{\mathrm{i}(\delta-\delta_{1})}=r_{1}re^{-\mathrm{i}(\delta-\delta_{1})}$. Since $r_{1},r\neq 0$, we may choose $\delta=\delta_{1}$ for $\arg(c_{1})$. We only need to compute $e^{\mathrm{i}\delta}$, so $e^{\mathrm{i}\delta}=\tilde{f}_{11}/|\tilde{f}_{11}|$. Now we have satisfied the conditions (16)–(17) for construction of a reflector that maps $\tilde{f}_{1}$ to $c_{1}e_{1}$, or to $-c_{1}e_{1}$.

We aim to compute $d$ or $s$ accurately. First, $w^{\ast}\widetilde{J}_{k}w$ is needed, where $w=\tilde{f}_{1}\pm f_{1}$ (with the addition for $s$ and the subtraction for $d$). Since $\tilde{f}_{1}$ and $f_{1}$ satisfy (16)–(17),

[TABLE]

Since (16) should be valid, we have $\operatorname{sgn}(\tilde{f}_{1}^{\ast}\widetilde{J}_{k}\tilde{f}_{1})=\tilde{\jmath}_{11}$, and both terms in the previous relation have the same sign. Therefore, to avoid unnecessary cancellation, our choice is $w=s$, and $H(s)$ from (15) is then equal to

[TABLE]

The update of a column $\tilde{f_{j}}$, $j>1$, is performed as $H(s)\tilde{f}_{j}=\tilde{f}_{j}+\tau s(s^{\ast}\widetilde{J}_{k}\tilde{f}_{j})$, which involves computing the $\widetilde{J}_{k}$-dot product of $s$ and $\tilde{f}_{j}$, scaling it by $\tau$, and calling the ZAXPY BLAS 1 routine. We opted for a sequential ZAXPY version, but it can be argued that for the extremely long columns a parallel version would be better suited.

The column updates are mutually independent, and therefore can be performed in a parallel-do loop over $j$, with the work being well balanced among the threads.

In our application it is not required that the reflector generators ($s_{k}$ and $\tau_{k}$ in the step $k$) are preserved, but that is nevertheless done with time overhead close to none in a separate complex matrix (the vector $s_{k}$ is stored in its $k$th column, with the leading rows set to 0) and a real vector (the scalar $\tau_{k}$ is stored as its $k$th element).

The situation is more complicated if a pair of pivot columns is chosen. According to [25], the reduction can be performed by the block Householder matrix defined by these two columns. However, the computation of a block Householder reflector is a more difficult approach than the computation of a variant of the URV factorization with $U$ hyperbolic, not unitary. The proposed factorization is similar, but not equal to the hyperbolic (sometimes also called signed) URV factorization presented in [30].

Definition 5.2**.**

*Let $F\in\mathbb{C}^{m\times n}$ and $J\in\mathbb{Z}^{m\times m}$, $J=\operatorname{diag}(\pm 1)$, be given matrices. Let $J^{\prime}=P^{T}JP$ for any permutation matrix $P$. A factorization $F=URV$, where $U\in\mathbb{C}^{m\times m}$ is $J^{\prime}$-unitary, $V\in\mathbb{C}^{n\times n}$ is unitary, and $R=\left[\begin{smallmatrix}R_{0}\\ 0\end{smallmatrix}\right]\in\mathbb{C}^{m\times n}$, with $R_{0}\in\mathbb{C}^{n\times n}$ upper triangular, is called a $(J,I)$ URV factorization according to $J$. *

Note that the $(J,I)$ URV factorization is not unique.

The Grammian matrix $A_{12}$ of the pivot block $\widetilde{F}_{12}\mathrel{\mathop{:}}=[\tilde{f}_{1}\ \tilde{f}_{2}]$, i.e., $A_{12}\mathrel{\mathop{:}}=\widetilde{F}_{12}^{\ast}\widetilde{J}_{k}\widetilde{F}_{12}$ (see Fig. 1) is nonsingular. Since the pivot strategy has chosen this block for the pivot block, the off-diagonal elements are larger in magnitude than the diagonal elements of $A_{12}$. The matrix $A_{12}$ will be diagonalized by a single Jacobi rotation $R_{k}$,

[TABLE]

This shows that the columns $[\hat{f}_{1}\ \hat{f}_{2}]\mathrel{\mathop{:}}=\widetilde{F}_{12}R_{k}$ are mutually $\widetilde{J}_{k}$ orthogonal – their scalar product is zero, and the squares of their $\widetilde{J}_{k}$ norms are $d_{k}$ and $d_{k+1}$, respectively. Since these columns are mutually $\widetilde{J}_{k}$ orthogonal, the only way to handle them is applying two successive hyperbolic Householder transformations to reduce the pivot matrix $\widetilde{F}_{12}$ to an upper triangular matrix $\widehat{F}_{12}$.

The first hyperbolic transformation, by $H(s_{k})$, reduces the first column to a single element. The second hyperbolic transformation, by $H(s_{k+1})$, then reduces the elements, already transformed by $H(s_{k})$, of the shortened (the first row of $\widehat{F}_{12}$ is not changed anymore) second column. We should then multiply $\widehat{F}_{12}$ by $R_{k}^{\ast}$ to obtain the reduced matrix $F_{12}$ with, generally full, topmost $2\times 2$ block

[TABLE]

Finally, the step counter $k$ is incremented by $2$ (instead of $1$ for a single pivot), and the process continues again with the column pivoting, as in subsection 5.2.

When present, the multiplications from right, first by the Jacobi rotation $R_{k}$, and then by $R_{k}^{\ast}$, cancel each other (since $R_{k}R_{k}^{\ast}=I_{2}$), so the right matrix $V$ in this URV-like factorization is in fact identity.

The Jacobi rotations are applied by calling the parallel ZROT routine, but should the columns be short enough, a sequential version of the routine could be warranted.

After the final step of the reduction (for $k=n$), the new $F$ and $J$ are square matrices of order $n$. More precisely, $J$ is the leading part of $\widetilde{J}$ as left after all permutations due to the row pivoting, and $F$, unlike the standard compact representation of the LAPACK’s QR factorizations, has zeros below the diagonal block set explicitly.

5.4 Testing

In Table 3 the wall execution times of both the JQR and the TSQR are shown. It is evident that there is still room for the future JQR’s efficiency improvement, which might be achieved by blocking and delaying the columns updates.

5.4.1 Prepermuting of $\widetilde{G}$

Table 3 also contains the wall time of preparing $\widetilde{G}P_{2}$ in parallel. The fastest way to prepermute the columns of $\widetilde{G}$ is to copy them to another matrix of the same size, with the column $j$ going to $\pi_{2}(j)$, where $\pi_{2}$ denotes the permutation represented by $P_{2}$. Such copying occurs in a parallel-do loop over $j$.

6 Phase 3 – generalized hyperbolic SVD

Vjeran Hari in his PhD thesis [9] developed a method for solving the generalized eigenproblem, when at least one of the two matrices is positive definite. The method is based on the ideas from the PhD thesis of Katharina Zimmermann [31], and has been revisited recently in [10, 11].

Based on the Hari–Zimmermann algorithm for the generalized eigenproblem, in [21] a one-sided method for computing the real generalized SVD has been derived. The main trick, how to obtain a one-sided method for the SVD from the two-sided method for the eigenproblem is always the same: think about the transformations in the two-sided fashion, and apply them from one (right or left) side on a matrix factor.

Since the elements of the pivot submatrices $\widehat{H}$ of $H$ and $\widehat{S}$ of $S$ are the scalar products of the columns of $F$ and $G$, respectively, it is easier to write the transformations in terms of the elements of $\widehat{H}$ and $\widehat{S}$,

[TABLE]

instead of in terms of the columns $f_{p}$, $f_{q}$, $g_{p}$, and $g_{q}$.

The original method consists of $3$ active transformations, and an auxiliary transformation that helps in coupling them all together.

6.1 Pointwise algorithm

The whole pointwise algorithm (with a pair of $2\times 2$ pivot submatrices in each annihilation step) is taken from the PhD thesis of Vjeran Hari [9]. However, we feel that the algorithm should be presented succinctly here, to aid its implementors, and also to incorporate some minor corrections.

6.1.1 Preprocessing

In the preprocessing step, $H$ and $S$ are scaled by a diagonal matrix $D$ such that $\operatorname{diag}(DSD)=I$, i.e.,

[TABLE]

Such preprocessing can be done only once, at the start of the algorithm, or it can be done before each annihilation step for the pivot column pair in question. The latter might seem redundant, but in a floating-point realization of the algorithm, after enough steps, $\operatorname{diag}(DSD)$ could veer off the identity enough to warrant such a rescaling. In that case, form a matrix $D_{0}$ that has, for the chosen pivot indices $(p,q)$, as its $p$th and $q$th diagonal entries $(s_{pp})^{-1/2}$ and $(s_{qq})^{-1/2}$, respectively, while being equal to the identity elsewhere. For the approach with a single prescaling, let $D_{0}=I$.

In both cases, let $\widehat{D}_{0}$ be a $2\times 2$ restriction of $D_{0}$ to the $p$th and the $q$th rows and columns. We have implemented the pivot pair prescaling in each annihilation step.

6.1.2 Diagonalization of $\widehat{S}_{0}$

In the first step the $2\times 2$ pivot submatrix $\widehat{S}_{0}$ of $S_{0}$ (at the crossings of the $p$th and the $q$th rows and columns) is diagonalized by a complex Jacobi rotation $\widehat{R}_{1}$, where $\widehat{R}_{k}$ for $k\geq 1$ is

[TABLE]

The same transformation is then applied to $H_{0}$, to keep the new pair equivalent to the original one. After that, the new pair is $(H_{1},S_{1})\mathrel{\mathop{:}}=(R_{1}^{\ast}H_{0}R_{1},R_{1}^{\ast}S_{0}R_{1})$, where $R_{1}=I$, except at the pivot positions, where $R_{1}=\widehat{R}_{1}$. Since $H$ and $S$ have been preprocessed as in subsection 6.1.1, the diagonal elements of $\widehat{S}_{0}$ are the same, and we may choose $\varphi_{1}=-\pi/4$ in (20). Now it is easy to determine that $\alpha_{1}=\arg(s_{pq})$. If $s_{pq}$, written in the trigonometric form, was $s_{pq}=xe^{\mathrm{i}\alpha_{1}}$, with $x=|s_{pq}|$, before the transformation, then after it we obtain

[TABLE]

6.1.3 Rescaling of $\widehat{S}_{1}$

The second step rescales the diagonal of $S_{1}$ to ones, and rescales $H_{1}$ with the same diagonal matrix. After the transformation, similar to the preprocessing step, we obtain $H_{2}\mathrel{\mathop{:}}=D_{2}H_{1}D_{2}$, $S_{2}\mathrel{\mathop{:}}=D_{2}S_{1}D_{2}$. From (21) we conclude that $\widehat{D}_{2}=\operatorname{diag}((1+x)^{-1/2},(1-x)^{-1/2})$. Elements of $\widehat{H}_{2}$ are

[TABLE]

where

[TABLE]

6.1.4 Diagonalization of $\widehat{H}_{2}$

In the third step the pivot submatrix $\widehat{H}_{2}$ of $H_{2}$ is diagonalized by a complex Jacobi rotation $\widehat{R}_{3}$ of the form (20). The third transformation can be written as $H_{3}=R_{3}^{\ast}H_{2}R_{3}$, $S_{3}=R_{3}^{\ast}S_{2}R_{3}$, where $R_{3}=I$, except at the pivot positions, where $R_{3}=\widehat{R}_{3}$. Then $\varphi_{3}$ in (20) is written as $\varphi_{3}=\vartheta+\pi/4$ to express the transformations in terms of $\vartheta$. The relations for the angles of a rotation that annihilates the off-diagonal element of $\widehat{H}_{2}$ are given by

[TABLE]

where $h=h_{qq}-h_{pp}$. The requirement $\gamma\mathrel{\mathop{:}}=\alpha_{3}-\alpha_{1}\in\left(-\pi/2,\pi/2\right]$ yields $\sigma=\operatorname{sgn}(h)$. This choice of $\sigma$ as the sign of $h$ and the constraint $\vartheta\in\left(-\pi/4,\pi/4\right]$ ensures the convergence of the algorithm (see [9]).

Since $\tan(2\vartheta+\pi/2)=-\cot(2\vartheta)=-1/\tan(2\vartheta)$, from (22) and (24), we obtain

[TABLE]

After the first three steps, the pivot submatrix $\widehat{S}_{3}$ is still diagonal (in fact identity), i.e., $\widehat{S}_{3}=\widehat{R}^{\ast}\widehat{S}_{0}\widehat{R}=I$, where

[TABLE]

This constructively shows that $\widehat{R}$ diagonalizes the matrix pair $(\widehat{H}_{0},\widehat{S}_{0})$.

6.1.5 Forming $\widehat{R}$ and $\widehat{Z}^{\prime}$

Hari in [9, Theorem 2.2] has proved the form of the general nonsingular matrix that diagonalizes a $2\times 2$ Hermitian positive definite matrix. The intention of representing $\widehat{R}$ as

[TABLE]

is to simplify (26). Comparing the elements of (26) and (27) we obtain

[TABLE]

Since in (28) we need only $\sin\gamma$ and $\cos\gamma$, from (24) it then follows

[TABLE]

The fourth step only deals with a formal simplification of $\widehat{R}$, by introducing a transformation $\Phi_{4}$ such that $H_{4}=\Phi_{4}^{\ast}H_{3}\Phi_{4}$, $S_{4}=\Phi_{4}^{\ast}S_{3}\Phi_{4}$. The matrix $\Phi_{4}$ is a diagonal matrix equal to identity, except at pivot positions, where $\widehat{\Phi}_{4}=\operatorname{diag}(e^{-\mathrm{i}\sigma_{p}},e^{-\mathrm{i}\sigma_{q}})$. Obviously, if $\widehat{R}$ diagonalizes the pair $(\widehat{H}_{0},\widehat{S}_{0})$, then the transformation $\widehat{Z}=\widehat{R}\widehat{\Phi}_{4}$ will leave the final diagonal matrices intact. Then, let $\widehat{Z}^{\prime}=\widehat{D}_{0}\widehat{Z}$.

6.1.6 Exceptional cases

There can be a few exceptions in the computations of the elements of the matrix $\widehat{Z}$ which have to be accounted for in the algorithm.

If $h_{pq}=s_{pq}=0$, we set $\widehat{Z}=I$, since both pivot submatrices are already diagonal. We could still apply the scaling by $\widehat{D}_{0}$, but not count that as a transformation.

If $h_{pq}\neq 0$, but $s_{pq}=0$ (i.e., $x=0$), we set $\alpha_{1}=0$ and proceed as described above to determine $\widehat{Z}$ as an ordinary Jacobi rotation that diagonalizes $\widehat{H}_{0}$.

If $h=v=0$ in (29), i.e., when $\arg{s_{pq}}=\arg{h_{pq}}$ and $h_{pp}=h_{qq}$, it can be shown that $\widehat{R}_{1}$, the Jacobi rotation that diagonalizes $\widehat{S}_{0}$, also diagonalizes $\widehat{H}_{0}$, so $\widehat{Z}=\widehat{R}_{1}\widehat{D}_{2}$.

6.1.7 Convergence criterion and finalization

As in [21], the convergence criterion in floating-point arithmetic has to take into account the relative magnitudes of the off-diagonal elements, compared to the diagonal ones, in the pivot pairs. Therefore, a pivot pair undergoes the transformation if, for the machine precision $\varepsilon$,

[TABLE]

If $\widehat{Z}^{\prime}$ turns out to be identity, it is not applied. If $\cos\varphi=\cos\psi=1$, such transformation is considered “small”, and “big” otherwise. Near the end of the process, the transformations turn out to be small, and (in a blocking variant, the last level of) the algorithm is stopped when no big transformations are encountered in a sweep, to avoid perpetually applying the transformations that spring only from the accumulation of the rounding errors. With blocking, the inner level(s) of the algorithm count all transformations (big and small) in a sweep for stopping. For details, see [18, 21].

Outputs

The algorithm stops when the columns of the in-place transformed $F^{\prime}$ and $G^{\prime}$ are numerically mutually $J$-orthogonal, and orthogonal, respectively. Let $Z^{\prime}$ be the accumulated product of the applied transformations. Then,

[TABLE]

and $U=F^{\prime}\Sigma_{F}^{\prime-1}$, $V=G^{\prime}\Sigma_{G}^{\prime-1}$, $\Sigma_{j}=\left((\Sigma_{F}^{\prime})_{j}^{2}+(\Sigma_{G}^{\prime})_{j}^{2}\right)^{1/2}$, $\Sigma=\operatorname{diag}(\Sigma_{1},\ldots,\Sigma_{n})$, $\Sigma_{F}=\Sigma_{F}^{\prime}\Sigma^{-1}$, $\Sigma_{G}=\Sigma_{G}^{\prime}\Sigma^{-1}$, $Z=Z^{\prime}\Sigma^{-1}$. For the inner levels of blocking, only the matrix $Z$ (and therefore also $\Sigma_{F}^{\prime}$, $\Sigma_{G}^{\prime}$, and $\Sigma$, but not $U$ and $V$) is required.

Accumulating $Z^{\prime-1}$

Optionally, $Z^{\prime-1}$ could be obtained by accumulating transformations $\widehat{Z}^{\prime}$ from the right and their inverses $\widehat{Z}^{\prime-1}$ from the left. Then, $X=Z^{-1}=\Sigma$, so Phase 4 would not be needed for the full G(H)SVD. This has not been implemented, since our main concern is a solution of the generalized eigenproblem.

6.1.8 The G(H)SVD algorithm

After obtaining the matrix $\widehat{Z}^{\prime}$, the pointwise implicit Hari–Zimmermann G(H)SVD algorithm can be written similarly to [21, Algorithm 3.1]. In Algorithm 2 we take into account the signature matrix $J$, since $H=F^{\ast}JF$, but with $J=I$ it reduces to a generalized SVD method.

If the prescaling as described in subsection 6.1.1 is performed only once, then in Algorithm 2 only $g_{p}^{\ast}g_{q}$ is computed, since the diagonal elements of $\widehat{S}$ are assumed to be unity. Otherwise, by passing once through the columns $g_{p}$ and $g_{q}$ all three dot products can be formed. Similar holds for $f_{p}$, $f_{q}$, and the three $J$-dot products of $\widehat{H}$.

6.2 Vectorization

Many computational building blocks of Algorithm 2 provide both the challenges and the opportunities for the SIMD vectorization. The most obvious such primitives are the $J$-dot products (and norms), which are computed combining the approaches presented in subsections 5.1 (for the unstructured patterns of signs in $J$) and 4.5.1 for a compact, partitioned representation of $J$ with only a number $n_{+}$ of the leading positive (and therefore $n-n_{+}$ tailing negative) signs.

6.2.1 ZVROTM

The column updates by $\widehat{Z}^{\prime}$ cannot be realized by a single call to a BLAS or a LAPACK routine (e.g., ZROT), since there are two angles involved in a transformation, and the columns are meant to be transformed in-place (overwritten). For that purpose, a ZVROTM routine has been implemented as a vectorized loop, that resembles a simplified version of the BLAS routine DROTM, but with the complex sines.

6.2.2 Transformations

The greatest challenge lies in computing the transformations in the SIMD-parallel way, where each vector lane $i$ computes $\widehat{Z}_{i}^{\prime}$ for its own pivot pair with indices $(p_{i},q_{i})$, as in [19] for the rotations in the Jacobi SVD method. Assume that $\mathsf{v}$ pivot pairs $(\widehat{H}_{i},\widehat{S}_{i})$ have been obtained, with their corresponding $p_{i}$ and $q_{i}$ indices all different. If there are fewer than $\mathsf{v}$ such pairs, let $\widehat{H}_{i}=\widehat{S}_{i}=I$ for the missing indices $i$. A vector (e.g., $\widehat{S}_{12}^{[i]}=g_{p_{i}}^{\ast}g_{q_{i}}$) has to be kept in an array of length $\mathsf{v}$, properly aligned in memory, in which the $i$th position holds data for the $i$th lane. To help the compiler, the complex arithmetic operations have been written in terms of the real and the imaginary parts of the complex numbers in the vectorizable regions.

For the start, check for which lanes their pivot pair has to be transformed by evaluating the criterion (30) in each lane. To aid the compiler, (30) can be rewritten as a branch-free arithmetic expression that has a non-zero value if and only if the criterion is fulfilled. If those values constitute a zero vector, no transformation is required for any $i$, and a fresh set of pivot pairs (if any remain) should be considered.

Then, compute $\widehat{Z}_{i}^{\prime}$ unconditionally (i.e., for all $i$). The idea is that the unneeded computation comes at no cost, while its results can be discarded afterwards. However, the exceptional cases from subsection 6.1.6 should be carefully dealt with, since a naïve branching might spoil the vectorization opportunities. The logical conditions are therefore arithmetized, while halting on the arithmetic exceptions is suppressed.

We assume that the second argument of the intrinsics MIN and MAX is returned when their first argument is a NaN. Such a behavior is not mandated by the Fortran standard (contrary to fmin and fmax in C), but is checked for at runtime by our code.

Let $\text{{B(i)}}=\text{{(1-MAX(V(i)/V(i),0))*(1-MAX(H(i)/H(i),0))}}$ and note that B(i) is 1 if, in (29), $v_{i}=h_{i}=0$ (so $\tan\gamma_{i}=\mathtt{NaN}$), and also if $h_{i}=\pm\infty$ (due to an overflow) with $v_{i}=0$; otherwise, it is 0. The computation resumes regardless of the value of B(i). However, $\widehat{Z}_{i}^{\prime}$ in those lanes where $\text{{B(i)}}=1$ is useless, so a new, correct $\widehat{Z}_{i}^{\prime}$ is taken according to the rules of subsection 6.1.6 at the end, sequentially for all such $i$. This is the only situation when $\widehat{Z}_{i}^{\prime}$ is not computed in one go for all $i$, but it should occur rarely in practice. A branch-free exception handling is also requried when $\widehat{S}_{12}^{[i]}=0$, since in calculating the polar form \big{|}\widehat{S}_{12}^{[i]}\big{|}e^{\mathrm{i}\phi_{i}} of $\widehat{S}_{12}^{[i]}$ the divisions of the real and the imaginary parts by the absolute value result in NaNs. If $\cos\phi_{i}$ is denoted by C(i), then selecting a default of $\cos\phi_{i}=1$ can be done by presetting the variable to 1, and taking $\text{{C(i)}}=\text{{MIN(x(i),C(i))}}$, where x is obtained by the vector division, or the multiplication by the reciprocal. For $\sin\phi_{i}$, the default value can be also set to 1, and the correct value of 0 is obtained by multiplying the intermediate result by a variable, set in the process of checking (30), that is 0 if $\widehat{S}_{12}^{[i]}=0$, and 1 otherwise.

6.3 Parallelization

There are $N\mathrel{\mathop{:}}=n(n-1)/2$ pairs of indices $(p_{i},q_{i})$ such that they belong to a strictly upper triangle (i.e., $1\leq p_{i}<q_{i}\leq n$) of the square matrices of order $n$. Let $\mathsf{P}\mathrel{\mathop{:}}=\{(p_{i},q_{i}):1\leq i\leq N\}$ be a set of those index pairs. At most $\left\lfloor n/2\right\rfloor$ such pairs can be chosen from $\mathsf{P}$ so that all their indices are distinct.

Let $\mathsf{S}_{j}$, for some $j\geq 1$, be a set of at most $\left\lfloor n/2\right\rfloor$ index pairs, with all indices distinct. Then, the pivot pairs formed from the columns of $F$ and $G$ and indexed by the elements of $\mathsf{S}_{j}$ can be transformed concurrently. We call $\mathsf{S}_{j}$ the $j$th parallel Jacobi step, and a sequence of steps $\mathsf{S}\mathrel{\mathop{:}}=(\mathsf{S}_{1},\mathsf{S}_{2},\ldots,\mathsf{S}_{\bar{n}})$ a parallel (quasi-)cyclic Jacobi strategy if $\bigcup_{j}\mathsf{S}_{j}=\mathsf{P}$, and under assumption that the sequence is repeated forever in principle (or, in practice, until the convergence criteria are met). A strategy is called cyclic if its steps are mutually disjoint; else, it is called quasi-cyclic. In a cycle (also called a sweep) all pivot pairs are accessed (at least once, but maybe more under a quasi-cyclic strategy), in $\bar{n}\geq n-1$ steps. We aim for the steps as large as possible.

6.3.1 Parallel strategies

Two classes of the parallel Jacobi strategies were under test: the modified modulus strategy (MM), described in, e.g., [21], and a generalization of the Mantharam–Eberlein strategy [17] (ME), described in [18]. The former is a quasi-cyclic strategy with $\bar{n}=n$, but easily generated on-the-fly as the computation progresses. The latter is cyclic, attains $\bar{n}=n-1$, has provided a faster execution than MM to the one-sided Jacobi SVD, with more accurate results, but is not readily available for all even $n$, and its convergence has not yet been proven.

6.3.2 Bordering

For an odd $n$, no more than $(n-1)/2$ pairs fit into a step, so $\bar{n}\geq n$. We therefore consider the strategies for even $n$ only, with the steps of size $n/2$, and when necessary border the matrices by appending a zero column and a zero row, except for the new element at position $(n+1,n+1)$, which is set to unity.

6.3.3 Vector-Parallel algorithm

Let a step $\mathsf{S}_{j}$, with $\mathsf{k}_{1}\mathrel{\mathop{:}}=n/2$ index pairs, be given. Then, partition $\mathsf{S}_{j}$ into $\mathsf{k}_{\mathsf{v}}\mathrel{\mathop{:}}=\left\lceil\mathsf{k}_{1}/\mathsf{v}\right\rceil$ disjoint subsets $\mathsf{V}_{k}$, where each subset has at most $\mathsf{v}$ pairs, i.e., $\mathsf{S}_{j}=(\mathsf{V}_{1},\mathsf{V}_{2},\ldots,\mathsf{V}_{\mathsf{k}_{\mathsf{v}}})$. For each subset, the requirements for the vectorized computation of $\widehat{Z}^{\prime}$ as described in subsection 6.2 are satisfied.

Now, let $\mathsf{t}\geq 1$ be a number of available (OpenMP) threads. Then, $\mathsf{S}_{j}$ is traversed with a parallel-do loop over the subsets, where a thread takes a chunk of subsets to be processed independently, while the subsets within a chunk are handled in sequence.

A sweep of such Vector-Parallel (VP) variant is shown in Algorithm 3.

Note that each thread has to have a private set of vector variables, which is most easily done by reserving a $\mathsf{v}\times\mathsf{t}$ rank-2 array for each variable and making the $l$th thread access the $l$th column of such an array, where $l$ is a thread’s unique number, $1\leq l\leq\mathsf{t}$.

6.4 Blocking

To better exploit the memory hierarchy by keeping data in the cache(s) longer, a multilevel blocking principle can be applied, with the Level 1 being the pointwise algorithm in its VP variant. In the next, second level of the algorithm the block columns of width $\mathsf{w}\geq 1$ take place of the single columns, and the $2\times 2$ pivot pairs are replaced by $(2\mathsf{w})\times(2\mathsf{w})$ block pivots. The same principle can be applied recursively further (e.g., see [18]), but we consider only the Level 2 algorithms here.

6.4.1 Block-column partitioning

There are two ways a matrix can be partitioned into block columns. The first one is to prescribe $\mathsf{w}$, and then split the matrix into at least $\left\lceil n/\mathsf{w}\right\rceil$ block columns of width at most $\mathsf{w}$, bearing in mind that the number of block columns has to be even, as explained in subsection 6.3.2, and reducing the maximal width accordingly. For $\mathsf{w}\geq 2$, all block columns can be made to contain either $\mathsf{w}$ or $\mathsf{w}-1$ (but no less) columns by redistributing their individual widths.

The second way, and the one we have chosen to implement, is to query at run-time a number $\mathsf{t}$ of threads to be used. A thread $l$ is to be assigned one pair of block columns with the block indices $(\mathsf{p}_{\mathsf{j}l},\mathsf{q}_{\mathsf{j}l})$ in the $\mathsf{j}$th block step, so the maximal width $\mathsf{w}$ is computed as $\left\lceil n/(2\mathsf{t})\right\rceil$, with exactly $2\mathsf{t}$ block columns. The block column widths $\mathsf{w}_{\mathsf{i}}$ are non-increasing across the whole partition and are equal to either $\mathsf{w}$ or $\mathsf{w}-1$.

The blocking overhead can dominate the actual computation time for the matrices small enough, so we assume that $n>2\mathsf{t}$, and suggest a pointwise algorithm otherwise.

Each thread allocates a contiguous storage (in its own NUMA region, but visible to the other threads) for $2\mathsf{w}$ columns (two block columns of the maximal width) of $F$, and similarly for $G$ and $Z$. The same amount of memory, and of the same shape, is additionally allocated for the columns of the “shadow” matrices $\mathsf{F}$, $\mathsf{G}$, and $\mathsf{Z}$. Let, for the $l$th thread in the $\mathsf{j}$th step of the $\mathsf{s}$th block sweep, the contents of that storage be named $\mathsf{X}_{\mathsf{j}l}^{[\mathsf{s}]}\mathrel{\mathop{:}}=[\mathsf{X}_{\mathsf{p}_{\mathsf{j}l}}^{[\mathsf{s}]}\ \mathsf{X}_{\mathsf{q}_{\mathsf{j}l}}^{[\mathsf{s}]}]$, with $\mathsf{X}\in\{F,G,Z,\mathsf{F},\mathsf{G},\mathsf{Z}\}$. If $\mathsf{w}_{\mathsf{p}_{\mathsf{j}l}}=\mathsf{w}-1$, the columns of $\mathsf{X}_{\mathsf{p}_{\mathsf{j}l}}$ are placed from the second column of $\mathsf{X}_{\mathsf{j}l}$, while the columns of $\mathsf{X}_{\mathsf{q}_{\mathsf{j}l}}$ are always placed from the column $\mathsf{w}+1$ of $\mathsf{X}_{\mathsf{j}l}$. That way all the columns of a block column pair are stored contiguously and can be viewed by the BLAS routines as a single matrix.

The storage for the block pivots and for the workspaces, along with a copy of $J$ and the strategy tables, is also preallocated per thread, in MCDRAM if possible. There are two strategy tables; one for the outer, Level 2 Jacobi strategy $\mathsf{S}^{[2]}$, and one for the inner, Level 1 strategy $\mathsf{S}^{[1]}$, which do not have to belong to the same class. The tables are fully initialized, with all the (block) steps, before the start of the iterations.

The initial data for $F_{1l}^{[1]}$ and $G_{1l}^{[1]}$ is loaded from the block columns $(\mathsf{p}_{1l},\mathsf{q}_{1l})\in\mathsf{S}_{1}^{[2]}$, $Z_{1l}^{[1]}$ is initialized to a corresponding part of $I_{n}$, while $\mathsf{F}_{1l}^{[1]}$, $\mathsf{G}_{1l}^{[1]}$, and $\mathsf{Z}_{1l}^{[1]}$ are undefined.

6.4.2 Processing the block pivots

In a step $\mathsf{j}$ (of a block sweep $\mathsf{s}$, which index we omit from the superscripts of the matrices for simplicity when it is implied by the context), the $l$th thread shortens its $J$ to $\widehat{J}_{\mathsf{j}l}$ and its block column pairs from $F_{\mathsf{j}l}$ to $\widehat{F}_{\mathsf{j}l}$ such that $\widehat{H}_{\mathsf{j}l}\mathrel{\mathop{:}}=\widehat{F}_{\mathsf{j}l}^{\ast}\widehat{J}_{\mathsf{j}l}\widehat{F}_{\mathsf{j}l}=F_{\mathsf{j}l}^{\ast}JF_{\mathsf{j}l}$, and $G_{\mathsf{j}l}$ to $\widehat{G}_{\mathsf{j}l}$ such that $\widehat{S}_{\mathsf{j}l}\mathrel{\mathop{:}}=\widehat{G}_{\mathsf{j}l}^{\ast}\widehat{G}_{\mathsf{j}l}=G_{\mathsf{j}l}^{\ast}G_{\mathsf{j}l}$, where $\widehat{F}_{\mathsf{j}l}$ and $\widehat{G}_{\mathsf{j}l}$ are both square, of order $2\mathsf{w}$ or $2(\mathsf{w}-1)$ (the order in between the two is necessarily odd, so the bordering from subsection 6.3.2 should then be applied).

There are two ways to shorten the block column pairs. The more efficient one is to form $\widehat{H}_{\mathsf{j}l}$ and $\widehat{S}_{\mathsf{j}l}$ explicitly. To do that, $F_{\mathsf{j}l}$ is copied to $\mathsf{F}_{\mathsf{j}l}$, and then the rows of $\mathsf{F}_{\mathsf{j}l}$ are scaled by $J$. Here, it is beneficial to have $J$ in the compact, run-length encoded form. Then, a single (parallel, or in our case, sequential) ZGEMM call computes $\widehat{H}_{\mathsf{j}l}=F_{\mathsf{j}l}^{\ast}(J\mathsf{F}_{\mathsf{j}l})$ and stores it temporarily into $\widehat{G}_{\mathsf{j}l}$. Then, $\widehat{F}_{\mathsf{j}l}^{\ast}$ and $\widehat{J}_{\mathsf{j}l}$ (partitioned to the positive and the negative sign blocks) are obtained by the Hermitian indefinite factorization with complete pivoting (see section 4), $\widehat{H}_{\mathsf{j}l}=\widehat{F}_{\mathsf{j}l}^{\ast}\widehat{J}_{\mathsf{j}l}\widehat{F}_{\mathsf{j}l}$. The factorization should reveal if $\widehat{H}_{\mathsf{j}l}$ is rank deficient, in which case the process stops (the computation could be retried with the (J)QR approach, as below, if the input data has been preserved). Finally, $\widehat{F}_{\mathsf{j}l}^{\ast}$ is copied, with the transposition and the complex conjugation applied, to $\widehat{F}_{\mathsf{j}l}$.

A simpler procedure is used for $\widehat{G}_{\mathsf{j}l}$. It suffices to compute $G_{\mathsf{j}l}^{\ast}G_{\mathsf{j}l}$ by a single (parallel, or in our case, sequential) ZHERK call, store $\widehat{S}_{\mathsf{j}l}$ temporarily to $\widehat{Z}_{\mathsf{j}l}$, and perform either the diagonally-pivoted Cholesky factorization, or even the Hermitian indefinite factorization with complete pivoting, to obtain $\widehat{S}_{\mathsf{j}l}=\widehat{G}_{\mathsf{j}l}^{\ast}\widehat{G}_{\mathsf{j}l}$. If $\widehat{S}_{\mathsf{j}l}$ is rank deficient or indefinite, the algorithm stops with an error message. A similar treatment has been implemented for $\widehat{H}_{\mathsf{j}l}$, if $J=I$. Otherwise, $\widehat{G}_{\mathsf{j}l}^{\ast}$ is copied, with the transposition and the complex conjugation applied, to $\widehat{G}_{\mathsf{j}l}$. A ZLASET call finally initializes $\widehat{Z}_{\mathsf{j}l}$ to $I$.

The (J)QR approach, similar to Phase 2, should be more accurate and even necessary when the conditions of the matrices $\widehat{H}_{\mathsf{j}l}$ and/or $\widehat{S}_{\mathsf{j}l}$ are so large that factorizing them after forming them explicitly might fail. Both $F_{\mathsf{j}l}$ and $G_{\mathsf{j}l}$ should then be copied to $\mathsf{F}_{\mathsf{j}l}$ and $\mathsf{G}_{\mathsf{j}l}$, respectively, and the JQR factorization on $\mathsf{F}_{\mathsf{j}l}$, followed by the column prepermutation and the column-pivoted QR factorization on $\mathsf{G}_{\mathsf{j}l}$, should be performed to obtain $\widehat{F}_{\mathsf{j}l}$ with $\widehat{J}_{\mathsf{j}l}$, and $\widehat{G}_{\mathsf{j}l}$, respectively. This has not been implemented, though.

The block pivots $\widehat{F}_{\mathsf{j}l}$ with $\widehat{J}_{\mathsf{j}l}$ and $\widehat{G}_{\mathsf{j}l}$ are handed over to a version of the Level 1 (VP) algorithm to be transformed. This, single-threaded VP version executes in the contexts of the already running threads. The Level 1 algorithm can either fully (implicitly) diagonalize $\widehat{H}_{\mathsf{j}l}$ and $\widehat{S}_{\mathsf{j}l}$, or at least iterate until a reasonably high number of the inner sweeps has been attained ($\mathsf{C}_{\max}=30$), in which case we talk about the Full Block (FB) variant; or pass over the block pivots only a prescribed number of times, e.g., once ($\mathsf{C}_{\max}=1$), in the Block-Oriented (BO) variant (for more details on both, see [21]). The former variant corresponds to a full two-sided annihilation of the off-diagonal of $\widehat{H}_{\mathsf{j}l}$ and $\widehat{S}_{\mathsf{j}l}$, while the latter implicitly reduces their off-diagonal norms.

Also, the BO variant exhibits similar execution times for every call of the Level 1 routine across all threads in a block step, while those can vary significantly, due to data, among the threads in the FB variant. On some platforms (e.g., the GPUs), that does not pose a huge problem [18], but on the CPUs it can cause delays on the synchronization primitives (OpenMP barriers) required between the block steps [21].

In any case, the transformations applied in the Level 1 are accumulated in $\widehat{Z}_{\mathsf{j}l}$, while the counters of all and of “big” transformations are added atomically to the respective counters (shared among the threads) for the current block sweep. With distributed memory, such counters can be updated by MPI_Allreduce collective calls.

6.4.3 Updating and exchanging the block columns

If no transformations have been applied in the Level 1, $F_{\mathsf{j}l}$, $G_{\mathsf{j}l}$, and $Z_{\mathsf{j}l}$ are then copied to $\mathsf{F}_{\mathsf{j}l}$, $\mathsf{G}_{\mathsf{j}l}$, and $\mathsf{Z}_{\mathsf{j}l}$, respectively. Otherwise, $\mathsf{F}_{\mathsf{j}l}=F_{\mathsf{j}l}\widehat{Z}_{\mathsf{j}l}$, $\mathsf{G}_{\mathsf{j}l}=G_{\mathsf{j}l}\widehat{Z}_{\mathsf{j}l}$, and $\mathsf{Z}_{\mathsf{j}l}=Z_{\mathsf{j}l}\widehat{Z}_{\mathsf{j}l}$ (three ZGEMM calls).

All computation have thus far been local to a thread, apart from the atomic operations. Now, by looking to $\mathsf{S}_{\mathsf{j}+1}^{[2]}$, it is easy to figure out which of the two block columns should be retained by the thread, and to which thread the other block column has to be sent (in fact, the whole communication pattern has been precomputed). Also, the block columns can swap roles: the first one in the current step can become the second one, in the presently owning or in the receiving thread, and vice versa, in the next step. A thread might even send away both its block columns (each to a different recipient), and receive two new block columns (each from a different sender).

It might seem unnecessary to perform the physical exchanges of data on the shared memory, but the reason behind them is twofold. First, most modern machines have their memory paritioned according to the speed of access by, or “proximity” to, each CPU (NUMA), and it is beneficial to bring the data close to (a thread bound to) the CPU that processes it. Second, such a design makes the algorithm convertible to a distributed memory one (e.g., by using an MPI process for what a thread does now).

To perform the block column copies, a thread has to know the memory addresses of the storage of the threads it communicates to. All such addresses are kept available to all threads, and the block columns from the “shadow” storage are copied by their present owners to the regular storage of their future ones. Before a block column pair can be copied (from or to), it has to be updated first, so there is an OpenMP barrier between the three updates above and the three copying actions (by two ZLACPY calls each, one per a block column). Also, a thread cannot continue with the next block step until it has the new data ready and its shadow storage available, which is enforced for all threads simultaneously by placing another barrier after the copying actions.

It is a legitimate question if it is worthwhile (and in what circumstances) to try hiding the communication behind the computation, maybe by relying on some tasking mechanism. For example, $\mathsf{F}_{\mathsf{p}_{\mathsf{j}l}}$ can be copied to $F_{\mathsf{q}_{\mathsf{j}+1,l^{\prime}}}$, and $\mathsf{F}_{\mathsf{q}_{\mathsf{j}l}}$ to $F_{\mathsf{p}_{\mathsf{j}+1,l^{\prime\prime}}}$, while $G_{\mathsf{j}l}$, $G_{\mathsf{j}l^{\prime}}$, and $G_{\mathsf{j}l^{\prime\prime}}$ are being updated. We leave those considerations for a future work.

Completing a sweep

At the end of a block sweep, the transformation counters are read and reset. If no “big” transformations have been applied in any of the threads in the block sweep, the process stops, and the outputs, now including $\Lambda_{F}$, $\Lambda_{G}$, and $\Lambda$, are generated, piecewise per thread, as described in subsection 6.1.7 and section 3.

6.5 Testing

The aim of testing was to establish what algorithm variant to recommend for practice, as well as should Phase 2 be employed, and if so, when.

6.5.1 Blocked vs. pointwise algorithms

In Table 4 it is shown that the Level 2 (BO) algorithm is several times faster than both the Level 1 (VP) and the Level 2 (FB) algorithms, and that blocking in general gives a significant advantage over the pointwise approach. Also, having more threads, and thus the smaller block pivot orders, benefits the FB algorithm, since the speedup with 64 versus 32 threads is there more than twofold. On the contrary, the speedup with twice more threads is less than twofold for the BO algorithm, since the formation of the block pivots in the BO variant takes a bigger portion of the overall time, compared to the Level 1 inner iterations. Nevertheless, the Level 2 (BO) algorithm is our choice for Phase 3.

6.5.2 Phase2 benefits

Phase 3 might have been run directly on the $\widetilde{F}$, $\widetilde{J}$, and $\widetilde{G}$. In Table 5 it is shown that the preprocessing of the tall-and-skinny inputs by Phase 2 into the square ones for Phase 3 is preferred, both time-wise and sweep-wise, but that advantage diminishes as the inputs approach a square-like form ($m\gtrsim n$).

6.5.3 Generalized eigenvalues

Please see the Figures S.2 and S.3 in the supplementary material, depicting the generalized eigenvalues $\Lambda(\text{A1})$–$\Lambda(\text{A4})$ and $\Lambda(\text{B1})$–$\Lambda(\text{B4})$, computed with 64 threads by the Level 2 (BO) Phase 3 GHSVD after shortening in Phase 2. The eigenvalues obtained with 32 threads differ a few ulps at most.

6.5.4 Comparison with ZHEGV(D)

In Table 6 the wall times for the explicit formation of $H$ (by the $\widetilde{J}$-scaling of $\widetilde{F}$ and the ZGEMM matrix multiplication) and $S$ (by the ZHERK matrix multiplication), as well as for the LAPACK’s generalized Hermitian eigensolvers ZHEGV and ZHEGVD on $(H,S)$, left in MCDRAM when possible, are shown, alongside the speedups of this approach (i.e., of forming of $H$ and $S$ and then calling either ZHEGV or ZHEGVD, with the eigenvectors also computed) vs. the Level 2 (BO) Phase 3 GHSVD on the inputs shortened by Phase 2, with 64 and 48 threads.

Please see Table S.2 in the supplementary material for the results with 32 and 24 threads.

The results suggest that both the CPU’s clock and the size of its low-level caches play a crucial role in performance of our approach for the problems large enough.

The main advantage of our approach is its ability to compute the generalized eigenproblem accurately when the LAPACK-based one can fail, not least due to the squaring of the condition number $\mathop{\kappa_{2}}(\widetilde{G})$ in the forming of $S=\widetilde{G}^{\ast}\widetilde{G}$, i.e.,

[TABLE]

that consequently leads to a failure in the Cholesky factorization, or to an unacceptable inaccuracy in the generalized eigenvalues. We believe that the demonstrated slowdown is a reasonable tradeoff for a reliable backup alternative in those difficult cases.

7 Phase 4 – optional computation of the full G(H)SVD

In Phase 4 the right generalized singular vector matrix $X=Z^{-1}$ is computed by the LU factorization with complete pivoting, $Z=P^{T}LUQ^{T}$, from the LAPACK routine ZGETC2, followed by solving a linear system $ZX=I_{n}$ for $X$, with $Z$ factored as above, by calling the sequential routine ZGESC2 in a parallel do loop. Each of $n$ loop iterations solves for one column of $X$, and the iteration space is divided among the same number of OpenMP threads, also available to ZGETC2, as used for the previous phases.

Should the G(H)SVD of the tall-and-skinny factors be required, $\widetilde{Z}=P_{2}^{T}Z$ and $\widetilde{X}=\widetilde{Z}^{-1}$ could also be computed in Phase 4, together with $\widetilde{U}=\widetilde{Q}_{F}U$ and $\widetilde{V}=\widetilde{Q}_{G}V$.

7.1 Relative errors in the full GHSVD

To measure the accuracy of the Phase 3 algorithm, we can look at the Frobenius norm of the error in the obtained GHSVD decomposition, relative to the Frobenius norm of the original matrix, i.e.,

[TABLE]

For a detailed description of computing the relative errors, see the supplementary material, subsection S.3.1, and for the full accuracy results, see Table S.3 therein.

The relative errors in the decomposition of $F$ range from $1.25\cdot 10^{-13}$ to $3.52\cdot 10^{-13}$, and from $1.12\cdot 10^{-13}$ to $7.22\cdot 10^{-13}$, for the datasets A and B, respectively. In the decomposition of $G$, the relative errors are slightly lower, ranging from $9.47\cdot 10^{-14}$ to $8.23\cdot 10^{-13}$, and from $7.98\cdot 10^{-14}$ to $4.73\cdot 10^{-13}$ for the datasets A and B, respectively. The relative errors are similar regardless the number of threads used in Phase 3.

From the range of errors it can be concluded both that our datasets are not highly ill-conditioned, and that in those cases the Phases 3 and 4 do not behave erratically. A rigorous stability analysis of the GHSVD method remains open for the future work.

7.2 Comparison with ZGGSVD3

The Phase 3 algorithm can also be used for the ordinary GSVD by setting $J=I$. It is therefore reasonable to compare it to the ZGGSVD3 LAPACK routine (from the parallel Intel MKL) serving the same purpose.

7.2.1 Dataset

A dataset C, comprising five Hermitian matrix pairs, has been generated by a call to the ZLATMS LAPACK testing routine for each matrix of the full bandwidth, with its pseudorandom eigenvalues uniformly distributed in $(0,1)$. The matrix orders are $k\cdot 1000$, $1\leq k\leq 5$, and the full GSVD is required.

7.2.2 Timing results

In Table 7 the speedup of Phase 3 followed by Phase 4, both with 64 threads, vs. ZGGSVD3 is shown. The speedup with 32 threads reaches a lower peak for $n=5000$, 140 times, as shown in the supplementary material, Table S.4.

From the tests it is evident that Phase 3 (with Phase 4 if needed) forms a very competitive, highly parallel algorithm for the ordinary GSVD as well, while being adaptable to a wide variety of the modern high performance computing hardware.

8 Future work

A (multi-)GPU version of the implicit Hari–Zimmermann algorithm for the ordinary GSVD [19, 20] shows the promising results and complements the CPU implementation described herein. Therefore, a GPU implementation of the GHSVD and Phase 1 might be a practical companion to the present research.

Acknowledgments

The authors would like thank the anonymous referees for their suggestions on improving the contents and the presentation of this manuscript.