# The distribution of the multiplicative index
of algebraic numbers over residue classes

Pieter Moree, Antonella Perucca, Pietro Sgobba  Max-Planck-Institut für
Mathematik, Vivatsgasse 7, D-53111 Bonn, Germany<EMAIL_ADDRESS>University of Luxembourg, Department of Mathematics, 6 Avenue de la Fonte,
L-4364 Esch-sur-Alzette, Luxembourg<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Let $K$ be a number field and $G$ a finitely generated torsion-free subgroup
of $K^{\times}$. Given a prime $\mathfrak{p}$ of $K$ we denote by
$\operatorname{ind}_{\mathfrak{p}}(G)$ the index of the subgroup
$(G\bmod\mathfrak{p})$ of the multiplicative group of the residue field at
$\mathfrak{p}$. Under the Generalized Riemann Hypothesis we determine the
natural density of primes of $K$ for which this index is in a prescribed set
$S$ and has prescribed Frobenius in a finite Galois extension $F$ of $K$. We
study in detail the natural density in case $S$ is an arithmetic progression,
in particular its positivity.

###### Key words and phrases:

Reductions of algebraic numbers, multiplicative index and order, primes in
arithmetic progression, natural density

###### 2020 Mathematics Subject Classification:

Primary: 11R45; Secondary: 11A07, 11R44

##  1\. Introduction

The distribution of the multiplicative index of an integer seems to have been
first studied by Pappalardi [12] in 1995. Under the Generalized Riemann
Hypothesis (GRH) he provided asymptotic formulae for $\sum_{p\leq
x}f(\operatorname{ind}_{p}(g))$, for $f$ satisfying fairly mild restrictions
(here and in the sequel we denote the rational primes by $p$). This line of
investigation was continued in 2012 by Felix and Murty [4] and later by Felix
for higher rank in [2]. Given a set of integers $S$ and a natural number $g$,
in [4] it was proven that

$\pi_{g,S}(x):=\left\lvert\\{p\leq x:\operatorname{ind}_{g}(p)\in
S\\}\right\rvert=c_{g,S}\operatorname{Li}(x)+O\Big{(}\frac{x}{(\log
x)^{2-\epsilon}}\Big{)},$ (1)

where $c_{g,S}$ is a constant defined by a series whose terms depend on the
set $S$, $\operatorname{Li}(x):=\int_{2}^{x}dt/\log t$ denotes the logarithmic
integral and $\epsilon>0$ is arbitrary. It is a difficult problem to determine
whether $c_{g,S}$ is positive or not, cf. Felix [3]. The special case where
$S$ is an arithmetic progression was already considered by Moree [8, Thm. 5]
in 2005. For example, he proved Theorem 8 below in case $G=\langle g\rangle$,
$F=K=\mathbb{Q}$.

In this paper we consider the behavior of $\pi_{g,S}(x)$, with $\mathbb{Q}$
replaced by a number field $K$ and $g$ by a finitely generated torsion-free
subgroup $G$ of $K^{\times}$. Instead of over rational primes, we sum now over
primes $\mathfrak{p}$ of norm $\leq x$. Under GRH we establish in Theorem 1,
see Section 2, an asymptotic similar to (1), but with a weaker error term
depending on the rank of $G$. In Section 3 we then restrict to the case where
$S$ consists of integers in an arithmetic progression $a\bmod d$. In Theorem 8
we show that in this case the natural density can be expressed as a linear
combination of at most $\varphi(d^{\prime})-1$ Artin-type constants, with
$d^{\prime}=d/(a,d)$. The positivity of the density is studied in Section 4,
the numerical evaluation of the Artin-type constants in Section 5. In the
final section we demonstrate our results by determining the density for two
examples and compare the outcome with an experimental approximation.

Likewise we can wonder about the above questions, but for the multiplicative
order, rather than the index. As far as the authors know, these were first
studied by Chinen and Murata [1] for $d=4$, and a little later by Moree by a
simpler method. Both Chinen and Murata, and independently Moree, went on to
independently write various further papers (he surveyed his results in [11]).
Under an appropriate generalization of the Riemann Hypothesis it turns out
that the natural density of primes $p\leq x$ such that the multiplicative
order of $g$ modulo $p$ is congruent to $a\bmod d$ exists. Denote it by
$\delta_{g}(a,d)$ and the associated counting function by $N_{g}(a,d)(x)$. The
proof of the existence of $\delta_{g}(a,d)$ by Moree is based on the identity

$N_{g}(a,d)(x)=\sum_{t=1}^{\infty}\left\lvert\\{p\leq
x\,:\,\operatorname{ind}_{p}(g)=t,\leavevmode\nobreak\ p\equiv 1+ta\bmod
dt\\}\right\rvert.$

The average density of elements of order congruent to $a\bmod d$ in a field of
prime characteristic also exists, but is a much simpler quantity, see Moree
[7]. It has very similar features to $\delta_{g}(a,d)$.

In the special case where $d$ divides $a$, we are just asking for the density
of primes $p$ such that $d$ divides the multiplicative order of $g$ modulo
$p$. This density is much easier to deal with and turns out to be a rational
number. This can be proven unconditionally, see for example [10, 17, 18].

Ziegler [19], using the approach of Moree, was the first to study the order in
arithmetic progression problem in the setting of number fields. His work was
generalized by Perucca and Sgobba in [14, 15], who obtained in particular
uniformity results for the distribution of the order. It is expected that,
likewise, some uniformity also holds for the distribution of the index into
suitably related congruence classes, however at the moment it is not clear how
to obtain such a result. For example, it does not follow from [14, Cor. 5.2],
in spite of the fact that congruence conditions on both the order and the size
of the multiplicative group lead to congruence conditions on the index. We
leave this as a research direction and as an open problem to the reader.

##  2\. The existence of the density of primes with prescribed index and
Frobenius

Let $K$ be a number field, and $F/K$ a finite Galois extension. Let $G$ be a
finitely generated and torsion-free subgroup of $K^{\times}$ having positive
rank $r$. Our goal is to determine the density of the set $P$ of primes
$\mathfrak{p}$ of $K$ (defined in the next theorem) with prescribed index and
Frobenius. The notation $F,K,G$ and $r$ will be maintained throughout. We also
set $K_{m,n}:=K(\zeta_{m},G^{1/n})$ for $m\mid n$, and similarly for
$F_{m,n}$. Further we make use of the following usual notation: $\zeta_{n}$
denotes an $n$-th primitive root of unity, $\mu$ the Möbius function, and
$\varphi$ Euler’s totient function. We write $\log^{a}x$ as shorthand for
$(\log x)^{a}$, and $(a,b)$ for $\gcd(a,b)$.

We recall that Landau’s prime ideal theorem states that

$\left\lvert\left\\{\mathfrak{p}:\operatorname{N}\mathfrak{p}\leq
x\right\\}\right\rvert=\operatorname{Li}(x)+O_{K}(xe^{-c_{K}{\sqrt{\log
x}}}),$ (2)

where $c_{K}>0$ is a constant depending on $K$.

######  Theorem 1.

(under GRH). Let $K$ be a number field, and let $G$ be a finitely generated
and torsion-free subgroup of $K^{\times}$ of positive rank $r$. Let $F/K$ be a
finite Galois extension, and let $C$ be a union of conjugacy classes in
$\operatorname{Gal}(F/K)$. Let $S$ be a non-empty set of positive integers.
Define

$P:=\\{\mathfrak{p}:\operatorname{ind}_{\mathfrak{p}}(G)\in
S,\,\operatorname{Frob}_{F/K}(\mathfrak{p})\in C\\},$

where $\mathfrak{p}$ ranges over the primes of $K$ unramified in $F$ and for
which $\operatorname{ind}_{\mathfrak{p}}(G)$ is well-defined. We let $P(x)$ be
the number of prime ideals in $P$ of norm $\leq x$. We have the asymptotic
estimate

$P(x)=\frac{x}{\log x}\sum_{t\in
S}\,\sum_{v=1}^{\infty}\frac{\mu(v)c(vt)}{[F_{vt,vt}:K]}+{O\left(\frac{x}{\log^{2-\frac{1}{r+1}}x}\right)},\\\
$

where

$c(n)=\left\lvert\\{\sigma\in\operatorname{Gal}(F_{n,n}/K):\sigma|_{K_{n,n}}=\operatorname{id},\,\sigma|_{F}\in
C\\}\right\rvert\,.$

The constant implied by the $O$-term depends only on $K$, $F$ and $G$.

This result in combination with the prime ideal theorem leads to the following
corollary.

######  Corollary 2.

(under GRH). Let $S$ be a non-empty set of positive integers. The natural
density of the primes $\mathfrak{p}$ of $K$ such that
$\operatorname{ind}_{\mathfrak{p}}(G)\in S$ and
$\operatorname{Frob}_{F/K}(\mathfrak{p})\in C$ exists and is given by

$\sum_{t\in S}\ \sum_{v=1}^{\infty}\frac{\mu(v)c(vt)}{[F_{vt,vt}:K]}\,.$

We will now formulate some preliminaries required for the proof of Theorem 1.
Our starting point is [14, Proposition 5.1], which was established for rank 1
in [19, Proposition 1].

######  Theorem 3.

(under GRH). For $x\geq t^{3}$, the number $R_{t}(x)$ of primes $\mathfrak{p}$
with norm up to $x$, unramified in $F$, and such that
$\operatorname{ind}_{\mathfrak{p}}(G)=t$ and
$\operatorname{Frob}_{F/K}(\mathfrak{p})\in C$ satisfies

$R_{t}(x)=\operatorname{Li}(x)\sum_{v=1}^{\infty}\frac{\mu(v)c(vt)}{[F_{vt,vt}:K]}+O\left(\frac{x}{\log^{2}x}\right)+O\left(\frac{x\log\log
x}{\varphi(t)\log^{2}x}\right).$

The constant implied by the $O$-term depends only on $K$, $F$ and $G$.

The following lemma is a straightforward generalisation of [8, Lemma 6],
taking into account that for every natural number $n$, the ratio

$C(n):=\frac{\varphi(n)n^{r}}{[K_{n,n}:K]}$ (3)

is bounded above by some constant $D$, depending only on $K$ and $G$ (see [14,
Theorem 1.1]).

######  Lemma 4.

For every real number $y\geq 1$ we have

$\sum_{t\leq
y}\,\sum_{n=1}^{\infty}\frac{\mu(n)}{[K_{nt,nt}:K]}=1+O\left(\frac{D}{y^{r}}\right),$

where the implied constant is absolute.

###### Proof.

We claim that

$\sum_{n>y}\frac{1}{n^{r}\varphi(n)}=O\Big{(}\frac{1}{y^{r}}\Big{)}.$

For $r=1$ this is due to Landau [6] who first proved that

$\sum_{n\leq x}\frac{1}{\varphi(n)}=A\log x+B+O\Big{(}\frac{\log
x}{x}\Big{)},$ (4)

with $A$ and $B$ explicit constants, and then applied partial integration. The
proof for arbitrary $r$ is completely analogous. Since
$\varphi(nt)\geq\varphi(n)\varphi(t)$, we obtain

$\sum_{t>y}\sum_{n=1}^{\infty}\frac{1}{[K_{nt,nt}:K]}\leq
D\sum_{t>y}\frac{1}{t^{r}\varphi(t)}\sum_{n=1}^{\infty}\frac{1}{n^{r}\varphi(n)}\ll\sum_{t>y}\frac{D}{t^{r}\varphi(t)}\ll\frac{D}{y^{r}},$
(5)

where we used that the fourth sum is bounded above by a constant not depending
on $r$. The estimate (5) shows that the double sum in the statement of Lemma 4
is absolutely convergent for all $y$. Thus, we may rearrange the double sum as
follows:

$\displaystyle\sum_{t=1}^{\infty}\sum_{n=1}^{\infty}\frac{\mu(n)}{[K_{nt,nt}:K]}$
$\displaystyle=\sum_{m=1}^{\infty}\sum_{s|m}\frac{\mu(m/s)}{[K_{m,m}:K]}=\sum_{m=1}^{\infty}\sum_{d|m}\frac{\mu(d)}{[K_{m,m}:K]}$
$\displaystyle=\sum_{m=1}^{\infty}\frac{1}{[K_{m,m}:K]}\sum_{d|m}\mu(d)=\frac{1}{[K_{1,1}:K]}=1,$

completing the proof. ∎

The following is a generalisation of Ziegler [19, Lemma 13].

######  Lemma 5.

(under GRH). We have

$\left\lvert\left\\{\mathfrak{p}:\operatorname{N}\mathfrak{p}\leq
x,\,\operatorname{ind}_{\mathfrak{p}}(G)>(\log
x)^{\frac{1}{r+1}}\right\\}\right\rvert=O\left(\frac{x}{\log^{2-\frac{1}{r+1}}x}\right),$

where the primes $\mathfrak{p}$ of $K$ are restricted to those for which
$\operatorname{ind}_{\mathfrak{p}}(G)$ is well-defined. The constant implied
by the $O$-term depends only on $K$ and $G$.

###### Proof.

The number of primes with ramification index or residue class degree at least
$2$ is of order $O(\sum_{p\leq\sqrt{x}}1)=O\left(\sqrt{x}/\log x\right)$. We
make use of the functions $R_{t}(x)$ from Theorem 3 with $F=K$. For any real
number $y\geq 1$, let $\mathcal{E}_{y}(x)$ be the number of primes
$\mathfrak{p}$ with $\operatorname{N}\mathfrak{p}\leq x$ and such that
$\operatorname{ind}_{\mathfrak{p}}(G)>y$. Notice that

$\mathcal{E}_{y}(x)=\left\lvert\left\\{\mathfrak{p}:\operatorname{N}\mathfrak{p}\leq
x\right\\}\right\rvert-\sum_{t\leq y}R_{t}(x)+O\left(\frac{\sqrt{x}}{\log
x}\right)\,.$

Landau’s prime ideal theorem (2) implies the (much) weaker estimate

$\left\lvert\left\\{\mathfrak{p}:\operatorname{N}\mathfrak{p}\leq
x\right\\}\right\rvert=\operatorname{Li}(x)+O\left(\frac{x}{\log^{2}x}\right)=\frac{x}{\log
x}+O\left(\frac{x}{\log^{2}x}\right),$ (6)

which is all we need. By Theorem 3 and Lemma 4 we obtain

$\displaystyle\sum_{t\leq y}R_{t}(x)$
$\displaystyle=\operatorname{Li}(x)\sum_{t\leq
y}\sum_{n=1}^{\infty}\frac{\mu(n)}{[K_{nt,nt}:K]}+O\left(\frac{xy}{\log^{2}x}\right)$
$\displaystyle+O\left(\frac{x\log\log x}{\log^{2}x}\sum_{t\leq
y}\frac{1}{\varphi(t)}\right)$
$\displaystyle=\operatorname{Li}(x)+O\left(\frac{x}{y^{r}\log
x}\right)+O\left(\frac{xy}{\log^{2}x}\right)$
$\displaystyle+O\left(\frac{x\log\log x}{\log^{2}x}\sum_{t\leq
y}\frac{1}{\varphi(t)}\right)\,.$

On taking $y=(\log x)^{1/(r+1)}$ we now obtain on invoking (6) and (4), the
estimate

$\mathcal{E}_{y}(x)=O\left(\frac{x}{y^{r}\log
x}\right)+O\left(\frac{xy}{\log^{2}x}\right)+O\left(\frac{x(\log\log
x)^{2}}{\log^{2}x}\right)=O\left(\frac{x}{\log^{2-\frac{1}{r+1}}x}\right),$

completing the proof. ∎

###### Proof of Theorem 1.

Set $\rho=2-\frac{1}{r+1}$. Lemma 5 with $y=(\log x)^{\frac{1}{r+1}}$ yields

$P(x)=\sum_{\begin{subarray}{c}t\leq y\\\ t\in
S\end{subarray}}R_{t}(x)+O\left(\frac{x}{\log^{\rho}x}\right).$

Estimating the sum as in the proof of Lemma 5 we obtain

$P(x)=\operatorname{Li}(x)\sum_{\begin{subarray}{c}t\leq y\\\ t\in
S\end{subarray}}\sum_{v=1}^{\infty}\frac{\mu(v)c(vt)}{[F_{vt,vt}:K]}+O\left(\frac{x}{\log^{\rho}x}\right)\,.$

Now we focus on the main term. We have

$\bigg{\rvert}\sum_{t\in
S}\sum_{v=1}^{\infty}\frac{\mu(v)c(vt)}{[F_{vt,vt}:K]}-\sum_{\begin{subarray}{c}t\leq
y\\\ t\in
S\end{subarray}}\sum_{v=1}^{\infty}\frac{\mu(v)c(vt)}{[F_{vt,vt}:K]}\bigg{\rvert}\leq\sum_{t>y}\sum_{v=1}^{\infty}\frac{1}{[F_{vt,vt}:F]}$

By (5) the right-hand side is bounded by $\ll y^{-r}$. Using this estimate the
proof is easily completed. ∎

##  3\. The distribution of the index over residue classes

Let $a,d$ be integers with $d\geq 2$. We study the density of primes
$\mathfrak{p}$ of $K$ such that $\operatorname{ind}_{\mathfrak{p}}(G)\equiv
a\bmod d$. Under GRH, by Theorem 1 this density exists and equals

$\operatorname{dens}_{G}(a,d):=\sum_{t\equiv a\bmod d}\ \sum_{v\geq
1}\frac{\mu(v)}{[K_{vt,vt}:K]}\,.$ (7)

The goal of this section is to prove Theorem 8, which expresses
$\operatorname{dens}_{G}(a,d)$ as a finite sum of terms depending on Dirichlet
characters $\chi$ of modulus $d$. These terms involve Artin-type constants
$B_{\chi}(r)$ that can be evaluated with multi-precision using Theorem 15,
thus allowing one to evaluate $\operatorname{dens}_{G}(a,d)$ with multi-
precision.

We start by explaining our notation. Given an integer $n\geq 1$ we let $G_{n}$
be the group of characters defined on $(\mathbb{Z}/n\mathbb{Z})^{\times}$, so
that $G_{n}\cong(\mathbb{Z}/n\mathbb{Z})^{\times}$. For a Dirichlet character
$\chi$ we denote by $h_{\chi}$ the (Dirichlet) convolution $\mu*\chi$ of the
Möbius function $\mu$ with $\chi$, that is $\mu*\chi(n)=\sum_{d\mid
n}\mu(d)\chi(n/d)$. Recall that the Dirichlet convolution of two
multiplicative functions is again a multiplicative function.

Put $w=\gcd(a,d),a^{\prime}=a/w$ and $d^{\prime}=d/w$. The integers $t\equiv
a\bmod d$ are of the form $wt^{\prime}$ with $t^{\prime}\equiv a^{\prime}\bmod
d^{\prime}$. Thus we can rewrite $\operatorname{dens}_{G}(a,d)$ as

$\operatorname{dens}_{G}(a,d)=\sum_{t\equiv a^{\prime}\bmod d^{\prime}}\
\sum_{v\geq 1}\frac{\mu(v)}{[K_{vwt,vwt}:K]}\,.$

This expression on its turn can be rewritten as

$\displaystyle\operatorname{dens}_{G}(a,d)$ $\displaystyle=\sum_{t\equiv
a^{\prime}\bmod d^{\prime}}\ \sum_{\begin{subarray}{c}v_{1}\geq 1\\\ t\mid
v_{1}\end{subarray}}\frac{\mu(v_{1}/t)}{[K_{v_{1}w,v_{1}w}:K]}$
$\displaystyle=\sum_{v_{1}\geq 1}\ \sum_{\begin{subarray}{c}t\equiv
a^{\prime}\bmod d^{\prime}\\\ t\mid
v_{1}\end{subarray}}\frac{\mu(v_{1}/t)}{[K_{v_{1}w,v_{1}w}:K]}$
$\displaystyle=\sum_{v_{1}\geq
1}\left(\frac{1}{\varphi(d^{\prime})}\sum_{\chi\in
G_{d^{\prime}}}\overline{\chi(a^{\prime})}h_{\chi}(v_{1})\right)\frac{1}{[K_{v_{1}w,v_{1}w}:K]}$
$\displaystyle=\frac{1}{\varphi(d^{\prime})}\sum_{\chi\in
G_{d^{\prime}}}\overline{\chi(a^{\prime})}\sum_{v_{1}\geq
1}\frac{h_{\chi}(v_{1})}{[K_{v_{1}w,v_{1}w}:K]}\,.$ (8)

In the second step we used that the double series is absolutely convergent
(see the proof of Lemma 4). In the third step we used [8, Lemma 9], where
$\chi$ runs over the Dirichlet characters modulo $d^{\prime}$.

We now focus on the final sum in (3). Recall the definition (3) of $C(n)$. By
Perucca et al. [16, Theorem 1.1] there exists an integer $n_{0}$ (depending
only on $G$ and $K$) such that

$C(n)=C(\gcd(n,n_{0}))\,.$ (9)

One can easily show that for $m\mid n$, one has $C(m)\mid C(n)$, and hence
$n_{0}$ can be taken to be the minimal integer satisfying

$C(n_{0})=\max_{n\geq 1}\frac{\varphi(n)n^{r}}{[K_{n,n}:K]}\,.$

By (9) we have

$\frac{1}{[K_{n,n}:K]}=\frac{C(\gcd(n,n_{0}))}{\varphi(n)n^{r}}\,,$

and therefore

$\sum_{n\geq 1}\frac{1}{[K_{n,n}:K]}=\sum_{g\mid
n_{0}}\sum_{\begin{subarray}{c}n\geq 1\\\
(n,n_{0})=g\end{subarray}}\frac{C(g)}{\varphi(n)n^{r}}\,.$

In our case,

$\sum_{v\geq 1}\frac{h_{\chi}(v)}{[K_{vw,vw}:K]}=\sum_{g\mid
n_{0}}\sum_{\begin{subarray}{c}v\geq 1\\\
(vw,n_{0})=g\end{subarray}}\frac{C(g)h_{\chi}(v)}{\varphi(vw)v^{r}w^{r}}\,.$
(10)

If $\sum_{v\geq 1}f(v)$ is some absolute convergent series, we have

$\displaystyle\sum_{\begin{subarray}{c}v\geq 1\\\
(vw,n_{0})=g\end{subarray}}f(v)$ $\displaystyle=\sum_{\begin{subarray}{c}v\geq
1\\\ (\frac{vw}{g},\frac{n_{0}}{g})=1\end{subarray}}f(v)=\sum_{v\geq
1}f(v)\sum_{n\mid\frac{n_{0}}{g},\,n\mid\frac{vw}{g}}\mu(n)$
$\displaystyle=\sum_{n\mid\frac{n_{0}}{g}}\mu(n)\sum_{\begin{subarray}{c}v\geq
1\\\
n\mid\frac{vw}{g}\end{subarray}}f(v)=\sum_{n\mid\frac{n_{0}}{g}}\mu(n)\sum_{\begin{subarray}{c}v\geq
1\\\ \frac{gn}{(gn,w)}\mid v\end{subarray}}f(v),$

where we used that $n$ divides the integer $vw/g$ if and only if $gn/(gn,w)$
divides $v$. Thus, in particular,

$\sum_{\begin{subarray}{c}v\geq 1\\\
(vw,n_{0})=g\end{subarray}}\frac{C(g)h_{\chi}(v)}{\varphi(vw)v^{r}w^{r}}=\frac{C(g)}{w^{r}}\sum_{n\mid\frac{n_{0}}{g}}\mu(n)\sum_{\begin{subarray}{c}v\geq
1\\\ \frac{gn}{(gn,w)}\mid
v\end{subarray}}\frac{h_{\chi}(v)}{\varphi(vw)v^{r}}\,.$

Inserting the right hand side into (10) and inserting the resulting expression
into (3) yields

$\operatorname{dens}_{G}(a,d)=\frac{1}{\varphi(d^{\prime})}\sum_{\chi\in
G_{d^{\prime}}}\overline{\chi(a^{\prime})}\sum_{g\mid
n_{0}}\frac{C(g)}{w^{r}}\sum_{n\mid\frac{n_{0}}{g}}\mu(n)\sum_{\begin{subarray}{c}v\geq
1\\\ \frac{gn}{(gn,w)}\mid
v\end{subarray}}\frac{h_{\chi}(v)}{\varphi(vw)v^{r}}\,.$

Denoting

$C_{\chi}(N,w,r)=\sum_{v\geq 1\atop N\mid
v}\frac{h_{\chi}(v)}{\varphi(vw)v^{r}},$

we can write this as

$\operatorname{dens}_{G}(a,d)=\frac{1}{\varphi(d^{\prime})}\sum_{\chi\in
G_{d^{\prime}}}\overline{\chi(a^{\prime})}\sum_{g\mid
n_{0}}\frac{C(g)}{w^{r}}\sum_{n\mid\frac{n_{0}}{g}}\mu(n)\,C_{\chi}\Big{(}\frac{gn}{(gn,w)},w,r\Big{)}\,.$
(11)

Let $\kappa(n)=\prod_{p\mid n}p$ denote the squarefree kernel of $n.$ Recall
that $h_{\chi}=\mu*\chi$. The following result is a special case of [8, Lemma
10] and expresses $C_{\chi}(N,w,r)$ as an Euler product.

######  Lemma 6.

We have

$C_{\chi}(N,w,r)=c_{\chi}(N,w,r)B_{\chi}(r),$

where

$c_{\chi}(N,w,r)=\frac{h_{\chi}(N)\kappa(Nw)}{N^{r+1}w}\prod_{p\mid
N}\frac{p^{r+1}}{p^{r+2}-p^{r+1}-p+\chi(p)}\prod_{p\nmid N\atop p\mid
w}\frac{p^{r+1}-1}{p^{r+2}-p^{r+1}-p+\chi(p)},$

and

$B_{\chi}(r)=\prod_{p}\Big{(}1+\frac{p(\chi(p)-1)}{(p-1)(p^{r+1}-\chi(p))}\Big{)}\,,$
(12)

where $p$ runs over all rational prime numbers.

######  Corollary 7.

We have $C_{\chi}(1,1,r)=\sum_{v\geq
1}\frac{h_{\chi}(v)}{v\varphi(v)}=B_{\chi}(r)$.

###### Proof of Lemma 6.

We distinguish two cases:
a) The case where $h_{\chi}(N)=0.$
We have to verify that $C_{\chi}(N,w,r)=0.$ Since $h_{\chi}$ is multiplicative
and we have $h_{\chi}(p^{k})=\chi(p)^{k-1}(\chi(p)-1),$ it follows that if
$h_{\chi}(N)=0,$ then there is a prime divisor $p$ of $N$ with $\chi(p)=1.$
Hence, $h_{\chi}(v)=0$ for all $v$ that are divisible by $N$ and so
$C_{\chi}(N,w,r)=0.$
b) The case where $h_{\chi}(N)\neq 0.$
We rewrite $C_{\chi}(N,w,r)$ as

$C_{\chi}(N,w,r)=\frac{h_{\chi}(N)}{\varphi(Nw)N^{r}}\sum_{v\geq
1}\frac{h_{\chi}(Nv)\varphi(Nw)}{h_{\chi}(N)\varphi(Nvw)v^{r}},$ (13)

and note that the argument is a multiplicative function in $v.$ We apply the
Euler product identity to evaluate the sum and obtain

$\prod_{p\mid N}\frac{p^{r+1}}{p^{r+1}-\chi(p)}\prod_{p\nmid N\atop p\mid
w}\Big{(}1+\frac{(\chi(p)-1)}{p^{r+1}-\chi(p)}\Big{)}\prod_{p\nmid
Nw}\Big{(}1+\frac{p(\chi(p)-1)}{(p-1)(p^{r+1}-\chi(p))}\Big{)},$

which can be rewritten as

$B_{\chi}(r)\prod_{p\mid
N}\frac{p^{r+1}(p-1)}{p^{r+2}-p^{r+1}-p+\chi(p)}\prod_{p\nmid N\atop p\mid
w}\frac{(p-1)(p^{r+1}-1)}{p^{r+2}-p^{r+1}-p+\chi(p)}\,.$

On inserting this in (13) and noting that

$\varphi(\kappa(Nw))=\prod_{p\mid N}(p-1)\prod_{p\nmid N\atop p\mid
w}(p-1),\quad\frac{\varphi(\kappa(Nw))}{\varphi(Nw)}=\frac{\kappa(Nw)}{Nw},$

the proof is completed. ∎

The densities we are interested in can be expressed as a finite linear
combination involving the constants $B_{\chi}(r)$. Our result generalizes [8,
Thm. 5] by Moree, who dealt with the case $F=K=\mathbb{Q}$ and $G$ of rank
$1$.

######  Theorem 8.

(under GRH). Let $a$ and $d$ be two natural numbers. Put $d^{\prime}=d/(a,d)$.
Assuming that the function $C(n)$, defined in (3), is explicitly given, we can
write

$\operatorname{dens}_{G}(a,d)=\sum_{\chi\in
G_{d^{\prime}}}d_{\chi}B_{\chi}(r),$

with the $d_{\chi}$ explicit complex numbers (they can be determined using
(11) and Lemma 6).

###### Proof.

In the identity (11) for $\operatorname{dens}_{G}(a,d)$ we make the
substitution

$C_{\chi}\bigg{(}\frac{gn}{(gn,w)},w,r\bigg{)}=c_{\chi}\bigg{(}\frac{gn}{(gn,w)},w,r\bigg{)}B_{\chi}(r)\,$

(which is allowed by Lemma 6). The constants $d_{\chi}$ are obtained by
factoring out the terms $B_{\chi}(r)$, so that for each $\chi\in
G_{d^{\prime}}$ we have

$d_{\chi}=\frac{\overline{\chi(a^{\prime})}}{\varphi(d^{\prime})}\sum_{g\mid
n_{0}}\frac{C(g)}{w^{r}}\sum_{n\mid\frac{n_{0}}{g}}\mu(n)c_{\chi}\Big{(}\frac{gn}{(gn,w)},w,r\Big{)}\,.\qed$

###  3.1. Generic aspects of the behaviour of $\operatorname{dens}_{G}(a,d)$

Generically the degree $[K_{vt,vt}:K]$ equals $vt\varphi(vt)$ if $G$ has rank
$1$. If every degree occurring in (7) would satisfy this, then we would obtain

$\rho(a,d):=\sum_{t\equiv a\bmod d}\ \sum_{v\geq
1}\frac{\mu(v)}{vt\varphi(vt)}\,.$

The inner sum is easily seen to equal $A\cdot r(t)$, with

$r(t)=\frac{1}{t^{2}}\prod_{p\mid t}\frac{p^{2}-1}{p^{2}-p-1}\,.$

Thus we can alternatively write

$\rho(a,d)=A_{1}\sum_{t\equiv a\bmod d}r(t)\,$

with

$A_{r}:=\prod_{p}\Big{(}1-\frac{1}{p^{r}(p-1)}\Big{)}$ (14)

the _rank $r$ Artin constant_. The “incomplete” rank $r$ Artin constant,
defined by restricting to $p$ odd, appears also in other works, such as in
Pappalardi [13]. For every $B>0$ we have, see [7, Theorem 4],

$\sum_{p\leq
x}\rho(p;a,d)=\rho(a,d)\operatorname{Li}(x)+O\Big{(}\frac{x}{\log^{B}x}\Big{)},$

with $\rho(p;a,d)$ the density of elements of in $\mathbb{F}_{p}^{*}$ having
index congruent to $a\bmod d$. Thus on average a finite field of prime order
has $\rho(a,d)$ elements having index congruent to $a\bmod d$. Two cases are
particularly easy.

######  Proposition 9 ([7, Proposition 4]).

One has

$\rho(0,d)=\frac{1}{d\varphi(d)}\quad\text{ and
}\quad\rho(d,2d)=\begin{cases}\rho(0,2d)&\text{ if }$d$\text{ is odd};\\\
3\rho(0,2d)&\text{ if }$d$\text{ is even}.\end{cases}$

In the remaining cases it is not difficult to express $\rho(a,d)$ in terms of
the $B_{\chi}(1)$’s, see [7, Prop. 6]. When $(a,d)=1$ this expression takes a
particularly simple form, namely

$\rho(a,d)=\frac{1}{\varphi(d)}\sum_{\chi\bmod
d}\overline{\chi(a)}B_{\chi}(1)\,.$ (15)

In the examples in Section 6 we will meet $\rho(a,d)$ again.

##  4\. The positivity of $\operatorname{dens}_{G}(a,d)$

As in the previous sections we consider a number field $K$, a finitely
generated and torsion-free subgroup $G$ of $K^{\times}$, and the natural
density $\operatorname{dens}_{G}(a,d)$ of the primes $\mathfrak{p}$ of $K$
such that $\operatorname{ind}_{\mathfrak{p}}(G)\equiv a\bmod d$. We are
interested in characterizing when this density is positive.

Recall that for a prime $\mathfrak{p}$ of $K$ of degree $1$ such that
$\operatorname{ind}_{\mathfrak{p}}(G)$ is well-defined, we have
$d\mid\operatorname{ind}_{\mathfrak{p}}(G)$ if and only if $\mathfrak{p}$
splits completely in $K_{d,d}$ (cf. [19, Lem. 2]). So by Chebotarev’s density
theorem we have

$\operatorname{dens}_{G}(0,d)=\frac{1}{[K_{d,d}:K]}>0,$ (16)

and thus we may suppose in the following that $0<a<d$.

We denote by $\operatorname{dens}_{G}(h)$ the density of primes $\mathfrak{p}$
such that $\operatorname{ind}_{\mathfrak{p}}(G)=h$ with $h$ a prescribed
integer, and by $n_{0}$ an integer satisfying $C(n)=C(\gcd(n,n_{0}))$, where
$C(n)$ was defined in (3). With this notation we are ready to recall the
following result by Järviniemi and Perucca:

######  Theorem 10 ([5, Main Thm. and Rem. 4.2], under GRH).

The density $\operatorname{dens}_{G}(h)$ is well-defined for all $h\geq 1$,
and we have $\operatorname{dens}_{G}(h)>0$ if and only if
$\operatorname{dens}_{G}(\gcd(h,n_{0}))>0$. For any set $S$ of positive
integers the following holds: if the density of primes $\mathfrak{p}$ of $K$
such that $\operatorname{ind}_{\mathfrak{p}}(G)\in S$ is positive, then there
is some $h\in S$ such that $\operatorname{dens}_{G}(h)>0$.

######  Proposition 11.

(under GRH). If $d\geq 2$ is coprime to $n_{0}$, then
$\operatorname{dens}_{G}(a,d)>0$.

###### Proof.

By Theorem 10 (taking $S$ to be the set of positive integers) we know that
there is some $h\geq 1$ such that $\operatorname{dens}_{G}(h)>0$. Moreover, we
deduce that there is an integer $h_{0}\mid n_{0}$ such that for every integer
$t$ coprime to $n_{0}$ we have $\operatorname{dens}_{G}(th_{0})>0$. We
conclude by taking $t\equiv 1\bmod n_{0}$ and $t\equiv ah_{0}^{-1}\bmod d$. ∎

The following result tells us in particular that for every prime number $\ell$
and for every $e\gg 0$ there is a positive density of primes $\mathfrak{p}$ of
$K$ such that $v_{\ell}(\operatorname{ind}_{\mathfrak{p}}(G))=e$.

######  Proposition 12.

For every prime number $\ell$ there is some non-negative integer $e_{\ell}$
(and we can take $e_{\ell}=0$ for all but finitely many $\ell$) such that for
every $e\geqslant e_{\ell}$ we have

$\operatorname{dens}_{G}(0,\ell^{e})>\operatorname{dens}_{G}(0,\ell^{e+1})\,.$

Under GRH, for every $n>0$ and for every integer $z$ we have
$\operatorname{dens}_{G}(z\ell^{e_{\ell}},\ell^{n})>0$.

###### Proof.

By Chebotarev’s density theorem the density of primes $\mathfrak{p}$ of $K$
such that $v_{\ell}(\operatorname{ind}_{\mathfrak{p}}(G))=e$ is given by
$1/[K_{\ell^{e},\ell^{e}}:K]-1/[K_{\ell^{e+1},\ell^{e+1}}:K]$, so the first
assertion follows from the eventual maximal growth of the Kummer degrees, see
[14, Lem. 3.2]. By the first assertion (and by applying Theorem 10 to the set
$S$ of positive integers having $\ell$-adic valuation equal to $e$) for every
$e\geq e_{\ell}$ there is some $b$ coprime to $\ell$ such that
$\operatorname{dens}_{G}(b\ell^{e})>0$. Then for every prime $q$ coprime to
$n_{0}$ we have $\operatorname{dens}_{G}(qb\ell^{e})>0$ so we may conclude by
selecting $q\equiv b^{-1}z\ell^{-v_{\ell}(z)}\bmod\ell^{n}$, which is possible
by Dirichlet’s theorem on primes in arithmetic progressions. ∎

If $x,y$ are positive integers, then we use the notation $\gcd(x,y^{\infty})$
to denote the positive integer obtained from $x$ by removing the prime factors
that do not divide $y$.

######  Theorem 13.

(under GRH). We have $\operatorname{dens}_{G}(a,d)>0$ if and only if

$\operatorname{dens}_{G}(a,\gcd(d,n_{0}^{\infty}))>0.$

###### Proof.

Set $d_{0}=\gcd(d,n_{0}^{\infty})$. The former inequality in the statement
clearly implies the latter because the integers congruent to $a\bmod d$ are
also congruent to $a\bmod d_{0}$. Now suppose that there is a positive density
of primes $\mathfrak{p}$ of $K$ such that
$\operatorname{ind}_{\mathfrak{p}}(G)\equiv a\bmod d_{0}$. From Theorem 10 we
deduce that there exists $h\geq 1$ such that $h\equiv a\bmod d_{0}$ and
$\operatorname{dens}_{G}(h)>0$. For every positive integers $t,s$ coprime to
$n_{0}$ such that $s\mid th$ we then have $\operatorname{dens}_{G}(th/s)>0$.
If we choose $t\equiv s\bmod d_{0}$, then $th/s\equiv a\bmod d_{0}$. We claim
that we may also choose $t,s$ so that $th/s\equiv a\bmod d/d_{0}$. Because of
the Chinese remainder theorem it will be possible to simultaneously ensure the
two conditions and hence $\operatorname{dens}_{G}(th/s)>0$ implies
$\operatorname{dens}_{G}(a,d)>0$. To prove the claim, we first choose $t,s$ so
that $\gcd(a,d/d_{0})=\gcd(th/s,d/d_{0})$ and then multiply $t$ by an integer
invertible modulo $d/d_{0}$ to obtain the requested congruence. ∎

######  Theorem 14.

(under GRH). The following conditions are equivalent:

  1. (1)

the density $\operatorname{dens}_{G}(a,d)$ is positive;

  2. (2)

there is an integer $A$ such that $A\equiv a\bmod d$ and
$\operatorname{dens}_{G}(\gcd(A,n_{0}))$ is positive.

Notice that it suffices to consider $A$ modulo $\operatorname{lcm}(d,n_{0})$.

###### Proof.

Write $D:=\operatorname{lcm}(d,n_{0})$. By Theorem 10 the density
$\operatorname{dens}_{G}(a,d)$ is positive if and only if there is an integer
$A\equiv a\bmod d$ for which $\operatorname{dens}_{G}(A,D)>0$. The latter
holds if and only if there is an index $h$ such that $h\equiv A\bmod D$ and
$\operatorname{dens}_{G}(h)>0$. Since $n_{0}\mid D$, we have
$\gcd(h,n_{0})=\gcd(A,n_{0})$, and hence by Theorem 10 we have that
$\operatorname{dens}_{G}(h)$ is positive if and only if
$\operatorname{dens}_{G}(\gcd(A,n_{0}))$ is positive.

The final assertion follows from the fact that the properties in (2) only
depend on $A$ modulo $\operatorname{lcm}(d,n_{0})$. ∎

##  5\. The Artin-type constants $B_{\chi}(r)$

Let $r\geq 1$ be an integer. Recall the Euler product definition (12) of
$B_{\chi}(r)$. For $r=1$ this was introduced in [7, Sec. 6] and denoted by
$B_{\chi}$, along with a variant $A_{\chi}$, where $p$ is restricted to those
primes for which $\chi(p)\neq 0$. We have

$B_{\chi}(1)=A_{\chi}\prod_{p\mid d}\left(1-\frac{1}{p(p-1)}\right),$

where $d$ is the modulus of the character. Note that $A_{\chi}=1$ in case
$\chi$ is the principal character.

If $\chi_{0}$ is the principal character, then $B_{\chi_{0}}(r)$ is a rational
number. This leaves at most $\varphi(d^{\prime})-1$ linearly independent
Artin-type constants, with $d^{\prime}=d/(a,d)$. For example, in case
$d^{\prime}=3$ and $d^{\prime}=4$ only one Artin-type constant is involved.
They are real numbers. As an illustration we point out the result that the
average density of elements of multiplicative order $\pm 1\bmod 3$ equals
$\frac{5}{16}\pm\frac{3}{10}B_{\chi_{3}}(1)$, where $\chi_{3}$ is the non-
principal character modulo $3$ and
$B_{\chi_{3}}(1)=\frac{5}{6}A_{\chi_{3}}=0.1449809353580\ldots$, see [7].

Approximating the numerical value of $B_{\chi}(r)$ by computing partial Euler
products, gives a quite poor accuracy. The following result allows us to do
rather better and generalizes [7, Thm. 6] to arbitrary $r$. It involves
special values of Dirichlet L-series. Recall that for $\Re(s)>1$ and $\chi$ a
Dirichlet character, we have

$L(s,\chi)=\sum_{n=1}^{\infty}\frac{\chi(n)}{n^{s}}=\prod_{p}\Big{(}1-\frac{\chi(p)}{p^{s}}\Big{)}^{-1}.$

######  Theorem 15.

Let $p_{1}(=2),p_{2},\ldots$ denote the sequence of consecutive primes and
$\chi$ be any Dirichlet character. Put

$\Lambda_{r}=A_{r}L(r+1,\chi)L(r+2,\chi)L(r+3,\chi).$

Then

$B_{\chi}(r)=E_{r,n}\,\Lambda_{r}\,\prod_{k=1}^{n}\left(1+\frac{\chi(p_{k})}{p_{k}(p_{k}^{r+1}-p_{k}^{r}-1)}\right)\Big{(}1-\frac{\chi(p_{k})}{p_{k}^{r+2}}\Big{)}\Big{(}1-\frac{\chi(p_{k})}{p_{k}^{r+3}}\Big{)}$

with

$1-\frac{1}{p_{n+1}^{r+2}}\leq|E_{r,n}|\leq 1+\frac{1}{p_{n+1}^{r+2}},$

provided that $r=1$ and $p_{n+1}\geq 5$, or $r=2$ and $p_{n+1}\geq 3$.

###### Proof.

Recall the definition (14) of $A_{r}$. Noting that

$(1-yt^{r+1})\left(\frac{1+\frac{(y-1)t^{r+1}}{(1-yt^{r+1})(1-t)}}{1-\frac{t^{r+1}}{1-t}}\right)=1+\frac{yt^{r+2}}{1-t-t^{r+1}},$

we obtain

$B_{\chi}(r)=A_{r}L(r+1,\chi)\prod_{k=1}^{\infty}\left(1+\frac{\chi(p_{k})}{p_{k}(p_{k}^{r+1}-p_{k}^{r}-1)}\right),$
(17)

on setting $y=\chi(p_{k})$ and $t=\frac{1}{p_{k}}$. We rewrite the infinite
product as

$L(r+2,\chi)L(r+3,\chi)\prod_{k=1}^{\infty}\left(1+\frac{\chi(p_{k})}{p_{k}(p_{k}^{r+1}-p_{k}^{r}-1)}\right)\Big{(}1-\frac{\chi(p_{k})}{p_{k}^{r+2}}\Big{)}\Big{(}1-\frac{\chi(p_{k})}{p_{k}^{r+3}}\Big{)}$

in order to improve its convergence. Denoting the $k$-th term in the infinite
product by $P_{r,k}$, we see that (17) holds with $E_{r,n}=\prod_{k\geq
n+1}P_{r,k}$.

It remains to estimate the relative error $E_{r,n}$ (which in general is a
complex number). Multiplying out

$(1-t-t^{r+1}+yt^{r+2})(1-yt^{r+2})(1-yt^{r+3})/(1-t-t^{r+1})$

gives

$1+\frac{yt^{r+4}}{1-t-t^{r+1}}\Big{(}1+t^{r-1}+(1-y)t^{r}-yt^{r+2}-yt^{2r+2}+y^{2}t^{2r+3}\Big{)},$

and leads to the estimate

$|P_{r,k}|\leq 1+t^{r+3}G(t,r),$

with

$G_{r}(t)=\frac{t(1+t^{r-1}+2t^{r}+t^{r+2}+t^{2r+2}+t^{2r+3})}{1-t-t^{r+1}}$

and $t=\frac{1}{p_{k}}$. Note that $G_{r}(t)$ is increasing in $t$ and
decreasing in $r$ in the region $0<t<1$ and $r\geq 1$. Thus

$|P_{r,k}|\leq 1+p_{k}^{-r-3}G_{r}(p_{k}^{-1})\leq
1+p_{k}^{-r-3}G_{r}(p_{n+1}^{-1})\text{\leavevmode\nobreak\
for\leavevmode\nobreak\ every\leavevmode\nobreak\ }k\geq n+1.$

As $t$ tends to zero, $G_{r}(t)$ tends to zero, and so we can choose $n$ so
large that $G_{r}(p_{n+1}^{-1})\leq 1$. Now

$|E_{r,n}|=\prod_{k\geq
n+1}|P_{r,k}|<\prod_{p>p_{n}}\left(1+\frac{1}{p^{r+3}}\right)<1+\sum_{m>p_{n}}\frac{1}{m^{r+3}}.$

Comparing the sum with an integral leads to the final estimate

$|E_{r,n}|\leq
1+\frac{1}{p_{n+1}^{r+3}}+\int_{p_{n+1}}^{\infty}\frac{dz}{z^{r+3}}\leq
1+\frac{1}{p_{n+1}^{r+2}},$

where the sum is over the integers $m>p_{n}$. Similarly,

$|E_{r,n}|>\prod_{p>p_{n}}\left(1-\frac{1}{p^{r+3}}\right)>1-\sum_{m>p_{n}}\frac{1}{m^{r+3}}>1-\frac{1}{p_{n+1}^{r+2}}.$

Some calculus shows that $G_{r}(\frac{1}{p})\leq 1$ if and only if $r=1$ and
$p\geq 5$ or $r\geq 2$ and $p\geq 3$. The proof is now completed on invoking
the if-part of this statement. ∎

######  Remark 16.

In the proof of [7, Thm. 6] there are a few typos:
For “$2+2t+t^{3}+t^{5}$” read “$2+2t+t^{3}+t^{4}+t^{5}$”.
For “$t\geq 127$” read “$t\leq 1/127$”.
For “$p_{n+!}$” read “$p_{n+1}$”.

##  6\. Two examples

In this section we demonstrate our results by two relatively easy, but
illustrative, examples for $K=\mathbb{Q}(\sqrt{5})$, $r=1$ and $d=5$. Some
examples for the same $r$ and $d$ values, but with $K=\mathbb{Q}$ are given in
Moree [9, Table 2].

$G$ | $\operatorname{dens}_{G}(0,5)$ | $\operatorname{dens}_{G}(1,5)$ | $\operatorname{dens}_{G}(2,5)$ | $\operatorname{dens}_{G}(3,5)$ | $\operatorname{dens}_{G}(4,5)$
---|---|---|---|---|---
$\frac{P_{0,5}(10^{6})}{\pi_{K}(10^{6})}$ | $\frac{P_{1,5}(10^{6})}{\pi_{K}(10^{6})}$ | $\frac{P_{2,5}(10^{6})}{\pi_{K}(10^{6})}$ | $\frac{P_{3,5}(10^{6})}{\pi_{K}(10^{6})}$ | $\frac{P_{4,5}(10^{6})}{\pi_{K}(10^{6})}$
$\langle\frac{1+\sqrt{5}}{2}\rangle$ | $0.100000$ | $0.418205$ | $0.296724$ | $0.0950872$ | $0.0899840$
$0.100093$ | $0.419351$ | $0.296954$ | $0.0947177$ | $0.0888838$
$\langle-\frac{5+\sqrt{5}}{2}\rangle$ | $0.100000$ | $0.451872$ | $0.266393$ | $0.0995570$ | $0.0821785$
$0.099787$ | $0.450979$ | $0.267518$ | $0.0996599$ | $0.0820564$
Table 1. Examples of densities $\operatorname{dens}_{G}(a,5)$ with
$K=\mathbb{Q}(\sqrt{5})$

In Table 1 we denote by $P_{a,d}(x)$ the number of primes $\mathfrak{p}$ of
$K$ of norm up to $x$ such that $\operatorname{ind}_{\mathfrak{p}}(G)\equiv
a\bmod d$, and by $\pi_{K}(x)$ the number of primes $\mathfrak{p}$ of $K$ with
norm up to $x$. The top row gives the theoretical density, the second row an
experimental approximation (both with rounding of the final decimal).

We will now treat these two examples without using the machinery of Section 3
(however, with complicated enough examples this becomes unavoidable). Our
approach requires some further notation. Given a divisor $\delta$ of an
integer $d_{1}$, we put

$\rho_{\delta,d_{1}}(a,d):=\sum_{t\equiv a\bmod d}\
\sum_{\begin{subarray}{c}v\geq 1\\\
(v,d_{1})=\delta\end{subarray}}\frac{\mu(v)}{vt\varphi(vt)}.$

###  6.1. First example

######  Proposition 17.

Set $K=\mathbb{Q}(\sqrt{5})$ and $G=\langle\frac{1+\sqrt{5}}{2}\rangle$. We
have

$\operatorname{dens}_{G}(0,5)=\frac{1}{10}=2\rho(0,5)$

and, for $1\leq a\leq 4$,

$\operatorname{dens}_{G}(a,5)=\frac{18}{19}\rho(a,5)=\frac{9}{38}\Big{(}\frac{19}{20}+\overline{\psi(a)}B_{\psi}(1)+\psi(a)B_{\psi^{3}}(1)+\psi^{2}(a)B_{\psi^{2}}(1)\Big{)}.$

###### Proof.

The first claim follows from (16) and Proposition 9. Next assume that $1\leq
a\leq 4$. We have $\rho(a,5)=\rho_{1,5}(a,5)+\rho_{5,5}(a,5)$. If $5\nmid t$,
then

$\sum_{5\mid
v}\frac{\mu(v)}{vt\varphi(vt)}=\sum_{w}\frac{\mu(5w)}{5wt\varphi(5wt)}=-\frac{1}{20}\sum_{5\nmid
w}\frac{\mu(w)}{wt\varphi(wt)}.$

We conclude that $\rho_{5,5}(a,5)=-\frac{1}{20}\rho_{1,5}(a,5)$. It thus
follows that $\rho_{1,5}(a,5)=\frac{20}{19}\rho(a,5)$ and
$\rho_{5,5}(a,5)=-\frac{1}{19}\rho(a,5)$. Since the degree $[K_{n,n}:K]$
equals $\varphi(n)n$ if $5\nmid n$ and $\frac{1}{2}n\varphi(n)$ otherwise, we
infer that
$\operatorname{dens}_{G}(a,5)=\rho_{1,5}(a,5)+2\rho_{5,5}(a,5)=\frac{18}{19}\rho(a,5)$.
The proof is completed on invoking (15) and noting that $\psi^{2}(a)$ is real
and $\overline{\psi^{3}(a)}=\psi(a)$. ∎

Approximations to $B_{\chi}(1)$ can be found in Table 2, where $\psi$ denotes
the character modulo $5$ determined uniquely by $\psi(2)=i$.

$\chi$ | $B_{\chi}(1)$
---|---
$\psi$ | $0.34645514515465\ldots+i\cdot 0.21283903970350\ldots$
$\psi^{2}$ | $0.12284254160167\ldots$
$\psi^{3}$ | $0.34645514515465\ldots-i\cdot 0.21283903970350\ldots$
$\psi^{4}$ | $0.95$
Table 2. The constants $B_{\chi}(1)$ for $d=5$

The character group has $\psi,\psi^{2},\psi^{3}$ and $\psi^{4}$ as elements,
with $\psi^{4}$ being the principal character. The table is taken from [7,
Table 3], where for $d\leq 12$ further approximations can be found. It was
kindly verified by Alessandro Languasco using Theorem 15 with $n=10^{6}$.

###  6.2. Second example

######  Proposition 18.

Set $K=\mathbb{Q}(\sqrt{5})$ and $G=\langle-\frac{5+\sqrt{5}}{2}\rangle$. Let
$1\leq a\leq 4$. One of $a$ and $a+5$ is even. Denoting this number by
$a_{1}$, we have

$\operatorname{dens}_{G}(a,5)=\frac{20}{19}\rho(a,5)-\frac{4}{19}\rho(a_{1},10).$

Furthermore, $\operatorname{dens}_{G}(0,5)=\frac{1}{10}$.

###### Proof.

Using (16) we see that $\operatorname{dens}_{G}(0,5)=\frac{1}{10}$. We will
determine $\operatorname{dens}_{G}(a,10)$ in case $5\nmid a$. The result then
follows on adding $\operatorname{dens}_{G}(a,10)$ and
$\operatorname{dens}_{G}(a+5,10)$.

Since
$\mathbb{Q}\Big{(}\sqrt{-\frac{5+\sqrt{5}}{2}}\Big{)}=\mathbb{Q}(\zeta_{5})$,
the degree $[K_{n,n}:K]$ equals $\varphi(n)n$ if $5\nmid n$, it equals
$\frac{1}{2}n\varphi(n)$ if $(n,10)=5$, and it equals $\frac{1}{4}n\varphi(n)$
if $10\mid n$. These degree considerations lead to

$\operatorname{dens}_{G}(a,10)=\begin{cases}\rho_{1,5}(a,10)+4\rho_{5,5}(a,10)&\text{
if }2\mid a;\\\ \rho_{1,5}(a,10)+2\rho_{5,10}(a,10)+4\rho_{10,10}(a,10)&\text{
if }2\nmid a.\end{cases}$

Reasoning as in the proof of Proposition 17 we deduce that
$\rho_{5,5}(a,10)=-\frac{1}{20}\rho_{1,5}(a,10)$ and
$\rho_{1,5}(a,10)=\frac{20}{19}\rho(a,10)$. It follows that
$\operatorname{dens}_{G}(a,10)=\frac{4}{5}\rho_{1,5}(a,10)=\frac{16}{19}\rho(a,10)$
in case $a$ is even.

If $a$ is odd, then so are the integers $t\equiv a\bmod 10$ and so
$\rho_{10,10}(a,10)=-\frac{1}{2}\rho_{5,10}(a,10)$, leading to
$\operatorname{dens}_{G}(a,10)=\rho_{1,5}(a,10)$. Reasoning as in the proof of
Proposition 17 we then deduce that
$\operatorname{dens}_{G}(a,10)=\frac{20}{19}\rho(a,10)$. ∎

For reasons of space we refrain here from explicitly writing out
$\operatorname{dens}_{G}(a,5)$ as a linear sum in the $B_{\chi}$’s, but we
will indicate how this is done. For $\rho(a,5)$ we use (15). For $a$ with
$5\nmid a$ we have by [7, Proposition 6] with $w=5$ and $\delta=2$,

$\rho(2a,10)=\frac{3}{8}\sum_{\chi\bmod
5}\overline{\chi(a)}\frac{B_{\chi}(1)}{2+\chi(2)}\,.$

## Acknowledgments

This project was started when the first author gave a talk in the Luxembourg
Number Theory Day 2019. He thanks the other authors for the invitation and for
several subsequent invitations. The second and third author thank the Max
Planck Institut für Mathematik and the first author for organizing a short
visit in October 2022. Thanks are also due to Alessandro Languasco for
verifying Table 2 using Theorem 15.

## References

  * [1] K. Chinen - L. Murata, _On a distribution property of the residual order of $a\,({\rm mod\leavevmode\nobreak\ }p)$_, I, II, J. Number Theory 105 (2004), 60–81, 82–100.
  * [2] T.A. Felix, _Higher rank generalizations of Fomenko’s conjecture_ , J. Number Theory 133 (2013), no. 5, 1738–1751.
  * [3] T.A. Felix, _The index of $a$ modulo $p$, SCHOLAR—a scientific celebration highlighting open lines of arithmetic research_, 83–96, Contemp. Math., 655, Centre Rech. Math. Proc., Amer. Math. Soc., Providence, RI, 2015.
  * [4] T.A. Felix - M. R. Murty, _A problem of Fomenko’s related to Artin’s conjecture_ , Int. J. Number Theory 8 (2012), no. 7, 1687–1723.
  * [5] O. Järviniemi - A. Perucca, _Unified treatment of Artin-type problems_ (2023), Res. Number Theory 9, no. 10, https://doi.org/10.1007/s40993-022-00418-6.
  * [6] E. Landau, _Über die zahlentheoretische Function $\varphi(n)$ und ihre Beziehung zum Goldbachschen Satz_, Göttinger Nachrichten, 177–186, 1900; Collected Works, Vol. I, pp. 106–115.
  * [7] P. Moree, _On the average number of elements in a finite field with order or index in a prescribed residue class_ , Finite Fields Appl. 10 (2004), no. 3, 438–463.
  * [8] P. Moree, _On the distribution of the order and index of $g\,({\rm mod\leavevmode\nobreak\ }p)$ over residue classes - I_, J. Number Theory 114 (2005), no. 2, 238–271.
  * [9] P. Moree, _On the distribution of the order and index of $g\,({\rm mod\leavevmode\nobreak\ }p)$ over residue classes - II_, J. Number Theory 117 (2006), no. 2, 330–354.
  * [10] P. Moree, _On primes $p$ for which $d$ divides $\operatorname{ord}_{p}(g)$_, Funct. Approx. Comment. Math. 33 (2005), 85–95.
  * [11] P. Moree, _On the distribution of the order over residue classes_ , Electron. Res. Announc. Amer. Math. Soc. 12 (2006), 121–128.
  * [12] F. Pappalardi, _On Hooley’s theorem with weights_ , Rend. Sem. Mat. Univ. Pol. Torino 53 (1995), no. 4, 375–388.
  * [13] F. Pappalardi, _On the $r$-rank Artin conjecture_, Math. Comp. 66 (1997), no. 218, 853–868.
  * [14] A. Perucca - P. Sgobba, _Kummer theory for number fields and the reductions of algebraic numbers_ , Int. J. Number Theory 15, no. 08 (2019), 1617–1633.
  * [15] A. Perucca - P. Sgobba, _Kummer theory for number fields and the reductions of algebraic numbers II_ , Uniform Distrib. Theory 15, no. 1 (2020), 75–92.
  * [16] A. Perucca - P. Sgobba - S. Tronto, _The degree of Kummer extensions of number fields_ , Int. J. Number Theory 17 (2021), no. 5, 1091–1110.
  * [17] K. Wiertelak, _On the density of some sets of primes - IV_ , Acta Arith. 43 (1984), no. 2, 177–190.
  * [18] K. Wiertelak, _On the density of some sets of primes $p$, for which $n|\operatorname{ord}_{p}(a)$_, Funct. Approx. Comment. Math. 28 (2000), 237–241.
  * [19] V. Ziegler, On the distribution of the order of number field elements modulo prime ideals, Uniform Distrib. Theory 1 (2006), no.1, 65–85.


# A new iterative method for solving a class of two-by-two block complex
linear systems

Davod Khojasteh Salkuyeh
Faculty of Mathematical Sciences, University of Guilan, Rasht, Iran
Center of Excellence for Mathematical Modelling, Optimization and
Combinational
Computing (MMOCC), University of Guilan, Rasht, Iran Corresponding author.
Emails<EMAIL_ADDRESS>

Abstract. We present a stationary iteration method, namely  Alternating
Symmetric positive definite and Scaled symmetric positive semidefinite
Splitting (ASSS), for solving the system of linear equations obtained by using
finite element discretization of a distributed optimal control problem
together with time-periodic parabolic equations. An upper bound for the
spectral radius of the iteration method is given which is always  less than 1.
So convergence of the ASSS iteration method is guaranteed. The induced ASSS
preconditioner is applied to accelerate the convergence speed of the GMRES
method for solving the system. Numerical results are presented to demonstrate
the effectiveness of both the ASSS iteration method and the ASSS
preconditioner.

Keywords: Iterative, PDE-constrained, optimization, convergence, finite
element, GMRES, preconditioning.
AMS Subject Classification: 49M25, 49K20, 65F10, 65F50.

##  1 Introduction

Consider the distributed control problems of the form: (see [14, 12, 24, 25])

$\displaystyle\min_{y,u}~{}~{}~{}\frac{1}{2}\int_{0}^{T}\int_{\Omega}|y(x,t)-y_{d}(x,t)|^{2}dxdt$
$\displaystyle+$
$\displaystyle\frac{\nu}{2}\int_{0}^{T}\int_{\Omega}|u(x,t)|^{2}dxdt,$
$\displaystyle\hskip
56.9055pt\text{s.t.}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{\partial}{\partial
t}y(x,t)-\Delta y(x,t)$ $\displaystyle=$ $\displaystyle u(x,t)~{}~{}{\rm
in}~{}~{}Q_{T},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\hskip
56.9055pt$ $\displaystyle y(x,t)$ $\displaystyle=$ $\displaystyle 0~{}~{}{\rm
on}~{}~{}\Sigma_{T},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\hskip
2.84544pt\hskip 56.9055pt$ $\displaystyle y(x,0)$ $\displaystyle=$
$\displaystyle y(x,T)~{}~{}{\rm
on}~{}~{}\partial\Omega,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\hskip
56.9055pt$ $\displaystyle u(x,0)$ $\displaystyle=$ $\displaystyle
u(x,T)~{}~{}{\rm
in}~{}~{}\Omega,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\hskip
1.42271pt\hskip 56.9055pt$

where $\Omega$ is an open and bounded domain in $\mathbb{R}^{d}$
($d\in\\{1,2,3\\})$ and its boundary $\partial\Omega$ is Lipschitz-continuous.
We introduce the space-time cylinder $Q_{T}=\Omega\times(0,T)$ and its lateral
surface  $\Sigma_{T}=\partial\Omega\times(0,T)$. Here, $\nu$ is a
regularization parameter, $y_{d}(x,t)$ is a desired state and $T>0$. Problems
of the form above arise widely in many areas of science and engineering [12,
13].

We may assume that $y_{d}(x,t)$ is time-harmonic, i.e.,
$y_{d}(x,t)=y_{d}(x)e^{i\omega t}$, with $\omega={2\pi k}/{T}$ for some
$k\in\mathbb{Z}$. If we substitute $y(x,t)$, $y_{d}(x,t)$ and $u(x,t)$ in the
problem then we get the following time-independent problem

$\displaystyle\min_{y,u}$
$\displaystyle\frac{1}{2}\int_{\Omega}|y(x)-y_{d}(x)|^{2}dx+\frac{\nu}{2}\int_{\Omega}|u(x)|^{2}dx,$
$\displaystyle s.t:$ $\displaystyle{\color[rgb]{0,0,0}i\omega y(x)}-\Delta
y(x)=u(x)~{}~{}\text{in}~{}~{}\Omega,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\hskip 54.06006pt{\color[rgb]{0,0,0}y(x)=0},\hskip
17.07182pt\text{on}~{}~{}\partial\Omega.$

We assume that $V_{h}\subset H_{0}^{1}(\Omega)$ is an $m$-dimensional vector
space spanned by the basis$\\{\varphi_{1},\varphi_{2},\ldots,\varphi_{m}\\}$.
The subspace $V_{h}$ is utilized for computing both the functions $y$ and $u$.
By applying the approach of discretize-then-optimization (see [16]), the above
problem can be described as the following form [13]

$\begin{array}[]{lll}\displaystyle\min_{y,u}&&\frac{1}{2}(\bar{y}-\bar{y}_{d})^{*}M(\bar{y}-\bar{y}_{d})+\frac{\nu}{2}\bar{u}^{*}M\bar{u},\\\
s.t:&&{\color[rgb]{0,0,0}iwM}\bar{y}+K\bar{y}=M\bar{u},\end{array}$ (1)

where the discretized negative Laplacian is interpreted by the matrix
$K=(\int_{\Omega}\nabla\varphi_{i}.\nabla\varphi_{j}dx)\in\mathbb{R}^{m\times
m}$ (the stiffness matrix) and
$M=(\int_{\Omega}\varphi_{i}\varphi_{j}dx)\in\mathbb{R}^{m\times m}$ denotes
the mass matrix. Here, the vectors $\bar{y}$, $\bar{y}_{d}$, and $\bar{u}$
denote the coefficients of the basis functions in $V_{h}$. We define the
Lagrangian functional for the discretized problem as

${\cal
L}(\bar{y},\bar{u},\bar{p})=\frac{1}{2}(\bar{y}-\bar{y}_{d})^{*}M(\bar{y}-\bar{y}_{d})+\frac{\nu}{2}\bar{u}^{*}M\bar{u}+{\bar{p}}^{*}({\color[rgb]{0,0,0}i\omega
M}\bar{y}+K\bar{y}-M\bar{u}),$

with $\bar{p}$ being the Lagrange multiplier of the constraint. By using the
Lagrange multiplier technique, we set $\nabla{\cal
L}(\bar{y},\bar{u},\bar{p})=0$, which is equivalent to

$\begin{pmatrix}M&0&K-i\omega M\\\ 0&\nu M&-M\\\ K+i\omega
M&-M&0\end{pmatrix}\begin{pmatrix}{\bar{y}}\\\ {\bar{u}}\\\
{\bar{p}}\end{pmatrix}=\begin{pmatrix}M{\bar{y}_{d}}\\\ 0\\\ 0\end{pmatrix}.$
(2)

It is known that, Eq. (2) gives the necessary and sufficient conditions for
the existence of a solution for the problem (1). It follows from the second
equation that $\bar{u}=\frac{1}{\nu}\bar{p}$ and substituting $\bar{u}$ in the
third equation, gives the following complex system of equations

$\left\\{\begin{array}[]{rl}M\bar{y}+(K-i{\omega}M)\bar{p}&=M{\bar{y}_{d}},\\\
(K+i{\omega}M)\bar{y}-\frac{1}{\nu}M\bar{p}&=0.\end{array}\right.$

The above system can be equivalently rewritten as

$A{x}=\begin{pmatrix}M&\sqrt{\nu}(K-i\omega M)\\\ \sqrt{\nu}(K+i\omega
M)&-M\end{pmatrix}\begin{pmatrix}{\bar{y}}\\\
{\bar{q}}\end{pmatrix}=\begin{pmatrix}{\hat{y}_{d}}\\\ 0\end{pmatrix}=b,$ (3)

where $\bar{q}=\bar{p}/\sqrt{\nu}$ and $\hat{y}_{d}=M{\bar{y}_{d}}$.

Since the systems of the form (3) are of very large size, it is important to
employ iterative methods in incorporated with suitable preconditioners for
solving these systems. In [13], Krendl et al. presented the real block
diagonal ($\mathcal{P}_{BD}$) and the alternative indefinite
($\mathcal{P}_{AI}$) preconditioners

$\displaystyle\mathcal{P}_{BD}=\begin{pmatrix}M+\sqrt{\nu}(K+\omega M)&0\\\
0&M+\sqrt{\nu}(K+\omega M)\end{pmatrix},$ (4)
$\displaystyle\mathcal{P}_{AI}=\begin{pmatrix}0&M+\sqrt{\nu}(K-i\omega M)\\\
M+\sqrt{\nu}(K+i\omega M)&0\end{pmatrix},$ (5)

for the system (3). Zheng et al. in [25] presented an improved version of the
block-diagonal preconditioner $\mathcal{P}_{BD}$. Zheng et al. in [24]
designed the block alternating splitting (BAS) iteration method for solving
the system (3) which can be written as

$\left\\{\begin{array}[]{rl}(\alpha V+H_{1})x^{(k+\frac{1}{2})}&=(\alpha
V-S_{1})x^{(k)}+\mathcal{P}_{1}b,{}{}{}{}{}{}\\\ (\alpha
V+H_{2})x^{(k+1)}&=(\alpha
V-S_{2})x^{(k+\frac{1}{2})}+\mathcal{P}_{2}b,{}{}{}\end{array}\right.$

where $\alpha>0$, $V$ is a symmetric positive definite matrix,

$H_{1}=\begin{pmatrix}M&0\\\ 0&M\end{pmatrix},\quad
H_{2}=\begin{pmatrix}\sqrt{\nu}K&0\\\ 0&\sqrt{\nu}K\end{pmatrix}$

and

$\displaystyle S_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\omega^{2}\nu}\begin{pmatrix}-{i\omega\nu}K&{\sqrt{\nu}}K\\\
-{\sqrt{\nu}}K&{i\omega\nu}K\end{pmatrix},\qquad
S_{2}=\begin{pmatrix}i\sqrt{\nu}\omega M&-M\\\ M&-i\sqrt{\nu}\omega
M\end{pmatrix},$ $\displaystyle\mathcal{P}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\omega^{2}\nu}\begin{pmatrix}-I&-i\omega\sqrt{\nu}I\\\
i\omega\sqrt{\nu}I&-I\end{pmatrix},\qquad\mathcal{P}_{2}=\frac{1}{1+\omega^{2}\nu}\begin{pmatrix}0&I\\\
I&0\end{pmatrix}.$

They proposed using the matrix $V=H_{1}$ as the preconditioner in practical
implementation. They proved that the BAS iteration method is convergent if
$\alpha\geq\nu\omega^{2}/2$. Numerical results presented in [24] demonstrate
that the BAS iteration method outperforms the GMRES method [19, 20]. On the
other hand, the BAS iteration method induces the preconditioner

$P_{BAS}(\alpha)=(1+\alpha)P(\alpha)\begin{pmatrix}\alpha M+\sqrt{\nu}K&0\\\
0&\alpha M+\sqrt{\nu}K\end{pmatrix},$

with

$P(\alpha)=\frac{1}{\alpha(2+\omega^{2}\nu)}\begin{pmatrix}I&(1+\omega^{2}\nu-i\omega\sqrt{\nu})I\\\
(1+\omega^{2}\nu+i\omega\sqrt{\nu})I&-I\end{pmatrix}.$

They also numerically illustrated that the parameter $\alpha=1+\nu\omega^{2}$
usually gives good results for the BAS iteration method and
$\alpha=(1+\nu\omega^{2})/(1+\sqrt{\nu}\omega)$ is a good choice for the BAS
preconditioner.  The reported numerical results in [24] showed that the BAS
preconditioner is superior to the preconditioners $\mathcal{P}_{BD}$ and
$\mathcal{P}_{AI}$.

Recently, Mirchi and Salkuyeh [15] stated the single block splitting (SBS)
iteration method for solving (3). They proved that if

$\frac{\nu\omega\mu_{\max}}{\sqrt{(1+\omega^{2}\nu)^{2}+\nu\mu_{\max}^{2}}}<1,$

then the SBS method is convergent, where $\mu_{\max}$ denotes the largest
eigenvalue of the matrix $Q=M^{-\frac{1}{2}}KM^{-\frac{1}{2}}$. However, if
$\sqrt{\nu}\omega$ is large, then the SBS method may converge slowly or even
fail to converge [15].

One may rewrite the system (3) in the real form

$\mathcal{K}\begin{pmatrix}r\\\ s\end{pmatrix}\equiv\begin{pmatrix}E&F^{T}\\\
F&-E\end{pmatrix}\begin{pmatrix}r\\\ s\end{pmatrix}=\begin{pmatrix}t\\\
0\end{pmatrix},$ (6)

where

$E=\begin{pmatrix}M&O\\\ O&M\end{pmatrix},\quad
F=\begin{pmatrix}\sqrt{\nu}K&-\sqrt{\nu}\omega M\\\ \sqrt{\nu}\omega
M&\sqrt{\nu}K\end{pmatrix},$ (7)

and

$r=\begin{pmatrix}\Re(\bar{y})\\\ \Im(\bar{y})\end{pmatrix},\quad
s=\begin{pmatrix}\Re(\bar{q})\\\ \Im(\bar{q})\end{pmatrix},\quad
t=\begin{pmatrix}\Re({\hat{y}_{d}})\\\ \Im({\hat{y}_{d}})\end{pmatrix}.$

It is noted that if $z$ is a complex vector, then $\Re(z)$ and $\Im(z)$ denote
the real and imaginary parts of $z$, respectively. Axelsson and Lukáš in [1]
(see also [4]) applied the preconditioned square block (PRESB) preconditioner

$\mathcal{C}=\begin{pmatrix}E+F+F^{T}&F^{T}\\\ F&-E\end{pmatrix},$ (8)

for system of linear equations of the form (6). They proved that every
eigenvalue of the matrix $\mathcal{C}^{-1}\mathcal{K}$ lies in the interval
$[\frac{1}{2},1]$ (see [1, 2, 3, 5]). Therefore, the GMRES method is indeed
appropriate for solving the preconditioned system. It turns out that for
applying the preconditioner $\mathcal{C}$, we need solving two systems where
the corresponding coefficient matrices are $E+F$ and $E+F^{T}$. We can also
solve these systems using the GMRES method incorporated with the PRESB
preconditioner. To apply the PRESB preconditioner for solving these systems,
we require to solve two subsystems with the coefficient matrix
$S=(1+w\sqrt{\nu})M+\sqrt{\nu}K$. Since the matrix $S$ is symmetric positive
definite, these systems can be solved either exactly by the Cholesky
factorization or inexactly by the conjugate gradient (CG) method. It is worth
nothing that if these systems are solved inexactly then we must employ the
flexible GMRES (FGMRES) [18] instead of the GMRES method.

In this paper, we present a stationary iteration method, namely Alternating
Symmetric positive definite and Scaled symmetric positive semidefinite
Splitting (ASSS), for solving the system (3) and verify its convergence
properties and the corresponding induced preconditioner.

We use the following notation in the rest of the paper. For a square matrix
$Z$, we denote $\sigma(Z)$ , $\rho(Z)$ and $\|Z\|_{2}$ for the spectrum,
spectral radius and 2-norm of the matrix, respectively. The Matlab notation
$(x;y)$ is used to denote the vector $(x^{T},y^{T})^{T}$. The imaginary unit
($\sqrt{-1}$) is shown by the letter $i$. The matrix $A\in\mathbb{R}^{n\times
n}$ is called positive real if $x^{T}Ax>0$, for every $0\neq
x\in\mathbb{R}^{n}$.

This paper is organized as follows. In Section 2 the ASSS iteration method is
introduced and its convergence analysis is given. Inexact version of the ASSS
iteration method is presented in Section 3. Section 4 is devoted to
introducing the ASSS preconditioner and its implementation issues. Numerical
results are presented in Section 5 to demonstrate the efficiency of the ASSS
iteration method and the corresponding preconditioner. Some concluding results
are drawn in Section 6.

##  2 The ASSS iteration method

Using the idea of [23], the system (3) can be written in the 4-by-4 block real
system

$\mathcal{A}{\bf
x}\equiv\left(\begin{array}[]{ccccc}M&0&\sqrt{\nu}K&\omega\sqrt{\nu}M\\\
0&M&-\omega\sqrt{\nu}M&\sqrt{\nu}K\\\ \sqrt{\nu}K&-\omega\sqrt{\nu}M&-M&0\\\
\omega\sqrt{\nu}M&\sqrt{\nu}K&0&-M\\\
\end{array}\right)\left(\begin{array}[]{l}\Re(\bar{y})\\\ \Im(\bar{y})\\\
\Re(\bar{q})\\\ \Im(\bar{q})\\\
\end{array}\right)=\left(\begin{array}[]{c}\Re(\hat{y}_{d})\\\
\Im(\hat{y}_{d})\\\ 0\\\ 0\end{array}\right)\equiv\hat{\mathbf{b}}.$ (9)

We define the matrices ${\mathcal{G}}_{1}$ and ${\mathcal{G}}_{2}$ as
following

${\mathcal{G}}_{1}=\left(\begin{array}[]{ccccc}I&0&0&\omega\sqrt{\nu}I\\\
0&I&-\omega\sqrt{\nu}I&0\\\ 0&-\omega\sqrt{\nu}I&-I&0\\\
\omega\sqrt{\nu}I&0&0&-I\\\
\end{array}\right),\quad{\mathcal{G}}_{2}=\sqrt{1+\nu\omega^{2}}\left(\begin{array}[]{ccccc}0&0&I&0\\\
0&0&0&I\\\ I&0&0&0\\\ 0&I&0&0\\\ \end{array}\right),\quad$

where $I\in\mathbb{R}^{m\times m}$ is the identity matrix. Then, the system
(9) can be equivalently rewritten as

$({\mathcal{G}}_{1}{\mathcal{M}}+\frac{\sqrt{\nu}}{\sqrt{1+\nu\omega^{2}}}{\mathcal{G}}_{2}\hat{{\mathcal{K}}}){\bf
x}=\hat{{\bf b}}.$ (10)

where

${\mathcal{M}}=\left(\begin{array}[]{cccc}M&0&0&0\\\ 0&M&0&0\\\ 0&0&M&0\\\
0&0&0&M\\\ \end{array}\right)\quad{\rm
and}\quad\hat{{\mathcal{K}}}=\left(\begin{array}[]{cccc}K&0&0&0\\\ 0&K&0&0\\\
0&0&K&0\\\ 0&0&0&K\\\ \end{array}\right).$

######  Lemma 1.

The matrix ${\mathcal{G}}_{1}$ is nonsingular and

${\mathcal{G}}_{1}^{-1}=\frac{1}{1+\nu\omega^{2}}\left(\begin{array}[]{ccccc}I&0&0&\omega\sqrt{\nu}I\\\
0&I&-\omega\sqrt{\nu}I&0\\\ 0&-\omega\sqrt{\nu}I&-I&0\\\
\omega\sqrt{\nu}I&0&0&-I\\\ \end{array}\right).$

Moreover, the matrix ${\mathcal{G}}:={\mathcal{G}}_{1}^{-1}{\mathcal{G}}_{2}$
is of the form

${\mathcal{G}}=\frac{1}{\sqrt{\nu(1+\nu\omega^{2})}}\left(\begin{array}[]{ccccc}0&\omega\nu
I&\sqrt{\nu}I&0\\\ -\omega\nu I&0&0&\sqrt{\nu}I\\\ -\sqrt{\nu}I&0&0&-\omega\nu
I\\\ 0&-\sqrt{\nu}I&\omega\nu I&0\\\ \end{array}\right).$

Furthermore, we have ${\mathcal{G}}^{2}=-I$ with $I$ being the identity matrix
of order $4m$, ${\mathcal{G}}^{-1}=-{\mathcal{G}}$ and
${\mathcal{G}}^{T}=-{\mathcal{G}}$.

###### Proof.

The proof is straightforward and is omitted here. ∎

Premultiplying both sides of the system (10) by ${\mathcal{G}}_{1}^{-1}$,
gives the system

$\mathcal{B}\bf{x}\equiv({\mathcal{M}}+{\mathcal{G}}{\mathcal{K}}){\bf x}={\bf
b},$ (11)

where $\mathcal{K}=\eta\hat{\mathcal{K}}$ with
$\eta=\sqrt{\nu}/\sqrt{1+\nu\omega^{2}}$ and ${\bf
b}={\mathcal{G}}_{1}^{-1}\hat{\bf{b}}$.

Given $\alpha>0$, we rewrite Eq. (11) as

$(\alpha I+{\mathcal{M}}){\bf x}=(\alpha I-{\mathcal{G}}{\mathcal{K}}){\bf
x}+{\bf b}.$ (12)

Eq. (11) can also be reformulated in the following form

$\mathcal{G}(\alpha I+{\mathcal{K}}){\bf
x}=(\alpha{\mathcal{G}}-\mathcal{M}){\bf x}+{\bf b}.$

Premultiplying both sides of this equation by $\mathcal{G}^{-1}$ and having in
mind that $\mathcal{G}^{-1}=-\mathcal{G}$, gives

$(\alpha I+{\mathcal{K}}){\bf x}=(\alpha I+\mathcal{G}\mathcal{M}){\bf
x}-\mathcal{G}{\bf b}.$ (13)

Now, using Eqs. (12) and (13) we establish the ASSS iteration method for
solving the system (11) as

$\left\\{\begin{array}[]{rl}(\alpha I+{\mathcal{M}}){\bf
x}^{(k+\frac{1}{2})}&=(\alpha I-{\mathcal{G}}{\mathcal{K}}){\bf x}^{(k)}+{\bf
b},\\\ (\alpha I+{\mathcal{K}}){\bf x}^{(k+1)}&=(\alpha
I+\mathcal{G}\mathcal{M}){\bf x}^{(k+\frac{1}{2})}-\mathcal{G}{\bf
b},\end{array}\right.$ (14)

where $\bf{x}^{(0)}$ is an initial guess.  It is worth noting that both of the
matrices ${\mathcal{M}}$ and ${\mathcal{K}}$ are symmetric positive definite
and the ASSS method is obtained using the splitting
$\mathcal{B}={\mathcal{M}}+{\mathcal{G}}{\mathcal{K}}$, which is indeed a
splitting of a symmetric positive definite matrix and a scaled (by
$\mathcal{G}$) symmetric positive definite matrix.

Eliminating ${\bf x}^{(k+\frac{1}{2})}$ from Eq. (14) yields the following
stationary iterative method

${\bf x}^{(k+1)}=\mathcal{T}_{\alpha}{\bf x}^{(k)}+\bf{f},$ (15)

where

$\mathcal{T}_{\alpha}=(\alpha I+\mathcal{K})^{-1}(\alpha
I+\mathcal{G}\mathcal{M})(\alpha I+\mathcal{M})^{-1}(\alpha
I-\mathcal{G}\mathcal{K}),$

is the iteration matrix of the ASSS method and

${\bf f}=\alpha(\alpha I+\mathcal{K})^{-1}(I-\mathcal{G})(\alpha
I+\mathcal{M})^{-1}{\bf b}.$

The following theorem states the convergence of the ASSS method.

######  Theorem 1.

Assume that the matrices $K$ and $M$ are symmetric positive definite, and let
$\alpha>0$. Then the spectral radius $\rho(\mathcal{T}_{\alpha})$ of the ASSS
iteration matrix

$\mathcal{T}_{\alpha}=(\alpha I+\mathcal{K})^{-1}(\alpha
I+\mathcal{G}\mathcal{M})(\alpha I+\mathcal{M})^{-1}(\alpha
I-\mathcal{G}\mathcal{K}),$

satisfies

$\rho(\mathcal{T}_{\alpha})\leq\gamma(\alpha)=\max_{\mu\in\sigma(M)}\frac{\sqrt{\alpha^{2}+\mu^{2}}}{\alpha+\mu}\max_{\lambda\in\sigma(K)}\frac{\sqrt{\alpha^{2}+(\eta\lambda)^{2}}}{\alpha+\eta\lambda},$

where $\eta=\sqrt{\nu}/\sqrt{1+\nu\omega^{2}}$. Hence, it holds that

$\rho(\mathcal{T}_{\alpha})\leq\gamma(\alpha)<1,\qquad\forall\alpha>0,$

which shows that the ASSS iteration method converges unconditionally.

###### Proof.

By setting $\hat{\mathcal{T}}_{\alpha}=(\alpha
I+\mathcal{K})\mathcal{T}_{\alpha}(\alpha I+\mathcal{K})^{-1}$, we see that

$\displaystyle\hat{\mathcal{T}}_{\alpha}$ $\displaystyle=$
$\displaystyle(\alpha I+\mathcal{G}\mathcal{M})(\alpha
I+\mathcal{M})^{-1}(\alpha I-\mathcal{G}\mathcal{K})(\alpha
I+\mathcal{K})^{-1}$ $\displaystyle=$
$\displaystyle\mathcal{R}_{\alpha}\mathcal{S}_{\alpha},$

where $\mathcal{R}_{\alpha}=(\alpha I+\mathcal{G}\mathcal{M})(\alpha
I+\mathcal{M})^{-1}$ and $\mathcal{S}_{\alpha}=(\alpha
I-\mathcal{G}\mathcal{K})(\alpha I+\mathcal{K})^{-1}$. Note that the matrices
$\hat{\mathcal{T}}_{\alpha}$ and ${\mathcal{T}}_{\alpha}$ are similar, and as
a result we deduce that

$\rho(\mathcal{T}_{\alpha})=\rho(\hat{\mathcal{T}}_{\alpha})=\rho(\mathcal{R}_{\alpha}\mathcal{S}_{\alpha})\leq\|\mathcal{R}_{\alpha}\mathcal{S}_{\alpha}\|_{2}\leq\|\mathcal{R}_{\alpha}\|_{2}\|\mathcal{S}_{\alpha}\|_{2}.$

On the other hand, it follows from $\mathcal{G}^{2}=-I$,
$\mathcal{G}^{T}=-\mathcal{G}$ and
$\mathcal{M}\mathcal{G}=\mathcal{G}\mathcal{M}$ that

$\displaystyle\|\mathcal{R}_{\alpha}\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\rho\left((\alpha I+\mathcal{M})^{-1}(\alpha
I-\mathcal{M}\mathcal{G})(\alpha I+\mathcal{G}\mathcal{M})(\alpha
I+\mathcal{M})^{-1}\right)$ (16) $\displaystyle=$ $\displaystyle\rho((\alpha
I+\mathcal{M})^{-1}(\alpha^{2}I+\alpha\mathcal{M}\mathcal{G}-\alpha\mathcal{G}\mathcal{M}-\mathcal{M}\mathcal{G}^{2}\mathcal{M})(\alpha
I+\mathcal{M})^{-1})$ $\displaystyle=$ $\displaystyle\rho\left((\alpha
I+\mathcal{M})^{-2}(\alpha^{2}I+\mathcal{M}^{2})(\alpha
I+\mathcal{M})^{-1}\right)$ $\displaystyle=$
$\displaystyle\max_{\mu\in\sigma(\mathcal{M})}\frac{\alpha^{2}+\mu^{2}}{(\alpha+\mu)^{2}}$
$\displaystyle=$
$\displaystyle\max_{\mu\in\sigma(M)}\frac{\alpha^{2}+\mu^{2}}{(\alpha+\mu)^{2}}.$

In the same way, we deduce that

$\|\mathcal{S}_{\alpha}\|_{2}^{2}=\max_{\lambda\in\sigma(K)}\frac{\alpha^{2}+(\eta\lambda)^{2}}{(\alpha+\eta\lambda)^{2}}.$
(17)

Since the matrices $M$ and $K$ are symmetric positive definite we have
$\mu,\lambda>0$, for all $\mu\in\sigma(M)$ and $\mu\in\sigma(K)$, and as
result we deduce that $\sqrt{\alpha^{2}+\mu^{2}}<\alpha+\mu$ and
$\sqrt{\alpha^{2}+(\eta\lambda)^{2}}<\alpha+\eta\lambda$. Therefore, we see
that $\|\mathcal{R}_{\alpha}\|_{2}<1$ and $\|\mathcal{S}_{\alpha}\|_{2}<1$.
Hence,
$\rho(\mathcal{T}_{\alpha})\leq\gamma({\alpha})=\|\mathcal{R}_{\alpha}\|_{2}\|\mathcal{S}_{\alpha}\|_{2}<1$,
which completes the proof. ∎

######  Remark 1.

Under the conditions of Theorem 1, we have

$\max_{\lambda\in\sigma(K)}\frac{\sqrt{\alpha^{2}+(\eta\lambda)^{2}}}{\alpha+\eta\lambda}<1,\quad\forall\alpha>0.$

Hence,

$\rho(\mathcal{T}_{\alpha})\leq\zeta(\alpha)=\max_{\mu\in\sigma(M)}\frac{\sqrt{\alpha^{2}+\mu^{2}}}{\alpha+\mu}.$
(18)

Similar to Corollary 2.1 in [6] the minimum value of the $\zeta({\alpha})$ is
obtained at

$\alpha^{*}=\sqrt{\mu_{\min}\mu_{\max}},$ (19)

where $\mu_{\min}$ and $\mu_{\max}$ are the smallest and largest eigenvalues
of the matrix $M$, respectively. It is worth noting that the parameter
$\alpha^{*}$ is independent of the parameters $\nu$ and $\omega$.  In
practice, we can use a few iterations of the power method for computing
$\mu_{\min}$ and the inverse power method to calculate $\mu_{\max}$.

######  Remark 2.

In Theorem 1 if the matrix $K$ is symmetric positive semidefinite, then

$\max_{\lambda\in\sigma(K)}\frac{\sqrt{\alpha^{2}+(\eta\lambda)^{2}}}{\alpha+\eta\lambda}\leq
1,\quad\forall\alpha>0,$

and we have $\rho(\mathcal{T}_{\alpha})\leq\zeta(\alpha)<1$, where
$\zeta(\alpha)$ was defined in Eq. (18). Hence, the ASSS iteration method is
still convergent unconditionally.

There are several methods to estimate the parameters of the iterative methods
and their induced preconditioners having similar structure to the ASSS
iteration method (see [11, 9, 22, 17, 8]). In [17], Ren and Cao proposed an
efficient estimation formula for the iteration parameter of a class of
alternating positive-semidefinite splitting preconditioner. The same idea was
used by Cao to estimate the parameter of the block positive-semidefinite
splitting preconditioner for generalized saddle point linear systems [8].
Similar to the strategy used in [17, 8], we can estimate the parameter of the
ASSS method. However, it is not as effective as the parameter
$\alpha^{*}=\sqrt{\mu_{\min}\mu_{\max}}$. We will shortly see in Section 5
that the value of $\alpha^{*}$ can be computed inexpensively.

##  3 Inexact version of ASSS

Let ${\bf x}^{(k)}$ be the computed solution solution at iteration $k$.
Setting

${\bf x}^{(k+\frac{1}{2})}={\bf x}^{(k)}+{\delta}^{(k)},$

and substituting it in the first step of Eq. (14) gives

$(\alpha I+{\mathcal{M}}){\delta}^{(k)}={\bf b}-{\cal B}{\bf x}^{(k)}=:{\bf
r}^{(k)}.$

In the same way if we set ${\bf x}^{(k+1)}={\bf
x}^{(k+\frac{1}{2})}+{\delta}^{(k+\frac{1}{2})}$ and substitute it in the
second step of Eq. (14), then we get

$\displaystyle{\color[rgb]{0,0,0}(\alpha
I+{\mathcal{K}}){\delta}^{(k+\frac{1}{2})}}$ $\displaystyle=$
$\displaystyle-{\mathcal{K}}{\bf
x}^{(k+\frac{1}{2})}+\mathcal{G}\mathcal{M}{\bf
x}^{(k+\frac{1}{2})}-\mathcal{G}{\bf b}$ $\displaystyle=$
$\displaystyle{\mathcal{G}}^{2}{\mathcal{K}}{\bf
x}^{(k+\frac{1}{2})}+\mathcal{G}\mathcal{M}{\bf
x}^{(k+\frac{1}{2})}-\mathcal{G}{\bf b}$ $\displaystyle=$
$\displaystyle-\mathcal{G}\left({\bf
b}-(\mathcal{M}+\mathcal{G}\mathcal{K}){\bf x}^{(k+\frac{1}{2})}\right)$
$\displaystyle=$ $\displaystyle-\mathcal{G}\left({\bf b}-\mathcal{B}{\bf
x}^{(k+\frac{1}{2})}\right)$ $\displaystyle=$ $\displaystyle-\mathcal{G}{\bf
r}^{(k+\frac{1}{2})}.$

Now using the above results we establish the inexact version of ASSS (IASSS)
as following.

Algorithm 1. The IASSS algorithm for solving $\mathcal{B}{\bf x}={\bf b}$.

  1. 1.

Choose an initial guess ${\bf x}^{(0)}$.

  2. 2.

For $k=0,1,2,\ldots$ until convergence, Do

  3. 3.

Compute ${\bf r}^{(k)}={\bf b}-\mathcal{B}{\bf x}^{(k)}$.

  4. 4.

Solve $(\alpha I+{\mathcal{M}}){\delta}^{(k)}={\bf r}^{(k)}$ approximately
using an iteration method.

  5. 5.

${\bf x}^{(k+\frac{1}{2})}:={\bf x}^{(k)}+\delta^{(k)}$.

  6. 6.

Compute ${\bf r}^{(k+\frac{1}{2})}={\bf b}-\mathcal{B}{\bf
x}^{(k+\frac{1}{2})}$.

  7. 7.

Solve $(\alpha I+\mathcal{K}){\delta}^{(k+\frac{1}{2})}=-\mathcal{G}{\bf
r}^{(k+\frac{1}{2})}$ approximately using an iteration method.

  8. 8.

${\bf x}^{(k+1)}:={\bf x}^{(k+\frac{1}{2})}+\delta^{(k+\frac{1}{2})}$.

  9. 9.

EndDo

In the steps 4 and 7 of Algorithm 1 two  linear systems of equations with the
coefficient matrices $\alpha I+{\mathcal{M}}$ and $\alpha I+{\mathcal{K}}$
should be solved. Since these matrices are SPD, they can be solved inexactly
using the CG method. However, since

$\alpha I+{\mathcal{M}}=\left(\begin{array}[]{cccc}\alpha I+M&0&0&0\\\
0&\alpha I+M&0&0\\\ 0&0&\alpha I+M&0\\\ 0&0&0&\alpha I+M\\\
\end{array}\right),$

we need to solve four systems with the same coefficient matrices $\alpha I+M$
and different right-hand sides. In fact if we set
${\hat{\delta}}^{(k)}=[{\delta}^{(k)}_{1};{\delta}^{(k)}_{2};{\delta}^{(k)}_{3};{\delta}^{(k)}_{4}]$
$\in\mathbb{R}^{4m}$ and ${\hat{\bf r}}^{(k)}=[{\bf r}^{(k)}_{1};{\bf
r}^{(k)}_{2};{\bf r}^{(k)}_{3};{\bf r}^{(k)}_{4}]\in\mathbb{R}^{4m}$ with
${\delta}^{(k)}_{i},{\bf r}^{(k)}_{i}\in\mathbb{R}^{m}$, for $i=1,2,3,4$, then
we only need to solve the linear system of equations with multiple right-hand
sides

$(\alpha I+M){\hat{\delta}}^{(k)}={\hat{\bf r}}^{(k)}.$

This system can be solved using the global CG algorithm with a suitable
preconditioner, which can be written as following. In this algorithm the inner
product used is $\langle X,Y\rangle=trace(X^{T}Y)$.

Algorithm 2. Global CG algorithm for $P^{-1}AX=P^{-1}B$.

  1. 1.

Compute $R_{0}=B-AX_{0}$, $Z_{0}=P^{-1}R_{0}$ and $P_{0}=Z_{0}$.

  2. 2.

For $j:=0,1,\ldots,$ until convergence, Do

  3. 3.

$\alpha_{j}:=\langle R_{j},Z_{j}\rangle/\langle AP_{j},P_{j}\rangle$.

  4. 4.

$X_{j+1}:=X_{j}+\alpha_{j}P_{j}$.

  5. 5.

$R_{j+1}:=R_{j}-\alpha_{j}AP_{j}$.

  6. 6.

$Z_{j+1}:=P^{-1}R_{j+1}$.

  7. 7.

$\beta_{j}:=\langle R_{j+1},Z_{j+1}\rangle/\langle R_{j},Z_{j}\rangle$.

  8. 8.

$P_{j+1}:=Z_{j+1}+\beta_{j}P_{j}$.

  9. 9.

Enddo

In the same way the system of Step 7 of Algorithm 1 can be solved. So in each
iteration of the IASSS iteration method two linear systems of equations with
multiple right-hand sides should be solved which can be solved using Algorithm
2.  The incomplete Cholesky factorization of the coefficient matrix can be
used as the preconditioner in Algorithm 2.

As Benzi and Golub mentioned in [7], adding $\alpha>0$ to the main diagonal of
$M$ (or $\eta K$) improves the condition number of the matrix. So this
improves, in turn, the convergence of the CG method much.

##  4 The ASSS preconditioner

It follows from $\mathcal{G}^{T}=-\mathcal{G}$ and
$\mathcal{G}\mathcal{K}=\mathcal{K}\mathcal{G}$ that

$\mathcal{B}^{T}+\mathcal{B}=\mathcal{M}+\mathcal{G}\mathcal{K}+\mathcal{M}+\mathcal{K}\mathcal{G}^{T}=\mathcal{M}+\mathcal{G}\mathcal{K}+\mathcal{M}-\mathcal{K}\mathcal{G}=2\mathcal{M},$

which shows that the matrix $\mathcal{B}$ is positive definite (positive
real). Hence, from Theorem 6.30 in [19] the restarted GMRES($\ell$)  for
solving the system (11) converges for any $\ell\geq 1$. On the other hand, if
we define

$\displaystyle\mathcal{P}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}(I+\mathcal{G})^{-1}(\alpha
I+\mathcal{M})\mathcal{G}(\alpha I+\mathcal{K}),$
$\displaystyle\mathcal{Q}_{\alpha}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha}(I+\mathcal{G})^{-1}(\alpha\mathcal{G}-\mathcal{M})(\alpha
I-\mathcal{G}\mathcal{K}),$

then

$\mathcal{B}=\mathcal{P}_{\alpha}-\mathcal{Q}_{\alpha}\quad\text{and}\quad\mathcal{T}_{\alpha}=\mathcal{P}^{-1}_{\alpha}\mathcal{Q}_{\alpha}.$

Hence, we get

$\mathcal{P}^{-1}_{\alpha}\mathcal{B}=I-\mathcal{T}_{\alpha},$

which shows that the eigenvalues of $\mathcal{P}^{-1}_{\alpha}\mathcal{B}$ are
clustered in a circle with radius 1 centered at $(1,0)$. Hence, the GMRES
method would be quite appropriate for solving the system incorporated with the
preconditioner $\mathcal{P}_{\alpha}$, i.e.,

$\mathcal{P}_{\alpha}^{-1}\mathcal{B}{\bf x}=\mathcal{P}_{\alpha}^{-1}{\bf
b}.$

In applying the preconditioner in each iteration of a Krylov subspace method
like GMRES a linear system of equations of the form $\mathcal{P}_{\alpha}s=r$
should be solved. Since

$\displaystyle\mathcal{P}_{\alpha}^{-1}$ $\displaystyle=$
$\displaystyle\alpha(\alpha I+\mathcal{K})^{-1}\mathcal{G}^{-1}(\alpha
I+\mathcal{M})^{-1}(I+\mathcal{G})$ $\displaystyle=$
$\displaystyle-\alpha(\alpha I+\mathcal{K})^{-1}\mathcal{G}(\alpha
I+\mathcal{M})^{-1}(I+\mathcal{G}),$

we can state the following algorithm for solving the system
$\mathcal{P}_{\alpha}s=r$.

Algorithm 3. Solution of $\mathcal{P}_{\alpha}s=r$.

  1. 1.

Compute $v=-\alpha(I+\mathcal{G})r$.

  2. 2.

Solve  $(\alpha I+\mathcal{M})w=v$ for $w$.

  3. 3.

Compute $z=\mathcal{G}w$.

  4. 4.

Solve  $(\alpha I+\mathcal{K})s=z$ for $s$.

Both of the systems in Steps 2 and 4 of the above algorithm can be solved
exactly using the Cholesky factorization of the matrices $\alpha I+M$ and
$\alpha I+\eta K$, or inexactly using the CG method. When these systems are
solved inexactly we can apply the global CG algorithm described in the
previous section. It is noted that, in this case the flexible version of GMRES
(FGMRES) [18] should be used for solving the main system instead of the GMRES
algorithm.

##  5 Numerical results

We consider the distributed control problem introduced in Section 1 in two-
dimensional case with the computational domain
$\Omega=(0,1)\times(0,1)\in\mathbb{R}^{2}$. The target state is set to be

$y_{d}(x,y)=\left\\{\begin{array}[]{cl}(2x-1)^{2}(2y-1)^{2},&\text{if}~{}~{}(x,y)\in(0,\frac{1}{2})\times(0,\frac{1}{2}),\\\
0,&\text{otherwise}.\end{array}\right.$ (20)

In our numerical test, we discretize the problem using the bilinear
quadrilateral Q1 finite elements with a uniform mesh [10]. Let
$D=\textrm{diag}(M)$, where $M$ is the mass matrix. In this case, the smallest
and largest eigenvalues of $N=D^{-1}M$ are $\frac{1}{4}$ and $\frac{9}{4}$,
respectively [21]. On the other hand, $D$ is a scalar multiplication of an
identity matrix, i.e., $D=\vartheta I$ for some $\vartheta>0$. Therefore, we
have $M=\vartheta N$. So, we deduce that the smallest and largest eigenvalues
of $M$ are $\mu_{\min}=\frac{1}{4}\vartheta$ and
$\mu_{\max}=\frac{9}{4}\vartheta$, respectively. Therefore, we have

$\alpha^{*}=\sqrt{\mu_{\min}\mu_{\max}}=\frac{3}{4}\vartheta.$

It is worth noting that the parameter $\alpha^{*}$ depends only on the mesh
size and independent of the parameters $\nu$ and $\omega$. In Table 1, we
disclose the values of $\vartheta$, $\mu_{\min}$, $\mu_{\max}$ and $\alpha$
for $h=2^{-k}$, $k=4,5,6,7$.

Table 1: The values of $\vartheta$, $\mu_{\min}$, $\mu_{\max}$ and $\alpha$
for $h=2^{-k}$, $k=4,5,6,7$.

$h$ | $\vartheta$ | $\mu_{\min}$ | $\mu_{\max}$ | $\alpha^{*}$
---|---|---|---|---
$2^{-4}$ | $1.7361\times 10^{-3}$ | $4.3403\times 10^{-4}$ | $3.9063\times 10^{-3}$ | $1.3021\times 10^{-3}$
$2^{-5}$ |  $4.3403\times 10^{-4}$ |  $1.0851\times 10^{-4}$ |  $9.7656\times 10^{-4}$ |  $3.2552\times 10^{-4}$
$2^{-6}$ | $1.0851\times 10^{-4}$ | $2.7127\times 10^{-5}$ | $2.4414\times 10^{-4}$ | $8.1380\times 10^{-5}$
$2^{-7}$ | $2.7127\times 10^{-5}$ | $6.7817\times 10^{-6}$ | $6.1035\times 10^{-5}$ | $2.0345\times 10^{-5}$

To generate the system (2) we have used the codes of the paper [16] which is
available at www.numerical.rl.ac.uk/people/rees/.  All runs are implemented in
Matlab R2017, equipped with a Laptop with 1.80 GHz central processing unit
(Intel(R) Core(TM) i7-4500), 6 GB RAM and Windows 7 operating system. For each
method, we report the number of iterations for the convergence and the elapsed
CPU time (in seconds). In the tables, a dagger (${\dagger}$) and a double
dagger (${\ddagger}$) mean that the method has not converged in 500 iterations
and 60 seconds, respectively.

We divide the numerical results in two parts. In the first part, we compare
the numerical results of the inexact version of the ASSS iteration method
described in Section 3 (denoted by IASSS) with those of the inexact version of
the BAS iteration method (denoted by IBAS). It is noted that in each iteration
of the BAS iteration method two subsystems with the coefficient matrix $\alpha
M+\sqrt{\nu}K$ and two systems with the coefficient matrix $M$ should be
solved. So in the two half-steps of the IBAS iteration method the subsystems
are solved using the global CG algorithm. The outer iteration is terminated as
soon as the residual norm of the system (3) is reduced by a factor of
$10^{6}$. The global CG iteration for solving the subsystems are stopped as
soon as the residual Frobenious norm of the residual matrix is reduced by a
factor of $10^{4}$. We always use a zero vector as an initial guess and the
maximum number of iterations is set to be 500.

In the IASSS iteration method we use the $\alpha^{*}$ computed by the formula
(19) (presented in Table 1), and in the IBAS iteration method the parameter
$\alpha$ is set to be $1+\nu\omega^{2}$ (as suggested in[24].) Numerical
results for different values of $h$, $\nu$ and $\omega$ have been presented in
Tables 2, 3 and 4. As we observe there is no significant difference between
the number iterations of the IASSS method when the parameters $h$, $\nu$ and
$\omega$ are changed. Comparing the numerical results of the IASSS method with
those of the IBAS method shows that the IBAS iteration method sometimes fails,
especially when $\nu$ and $\omega$ are large. However, this is not the case
for the IASSS iteration method. Nevertheless, when the values of $\omega$ and
$\nu$ are small enough, the CPU time of IBAS is often less than that of IASSS.

For the second part of our experiments, we compare numerically the performance
of the ASSS preconditioner (P-ASSS) with those of the BAS preconditioner
(P-BAS), the PRESB preconditioner (P-PRESB) and the preconditioner
$\mathcal{P}_{BD}$ (P-BD ). To do so, we use the flexible version of the GMRES
(FGMRES) method in conjunction with aforementioned preconditioners. In the
implementation of the PRESB preconditioner (8), the systems with the
coefficient matrices $E+F$ and $E+F^{T}$ are solved using the FGMRES method
incorporated with the PRESB preconditioner. The innermost subsystems with the
coefficient matrix $S=(1+w\sqrt{\nu})M+\sqrt{\nu}K$ in the PRESB
preconditioner and the systems with the coefficient matrix
$T=(1+\omega)M+\sqrt{\nu}K$ in $\mathcal{P}_{BD}$ are solved using the CG
method in conjunction with the incomplete Cholesky factorization with dropping
tolerance $0.001$ as a preconditioner. For applying the ASSS and the BAS
preconditioners all the subsystems are solved using the global CG method with
incomplete Cholesky factorization with dropping tolerance $0.001$ as a
preconditioner. We use $\alpha^{*}$ reported in Table 1 for the ASSS
preconditioner and $\alpha=(1+\nu\omega^{2})/(1+\sqrt{\nu}\omega)$ for the BAS
preconditioner (as suggested in [24]).

The iteration of the FGMRES method as the outer iteration is stopped as soon
as the residual norm is reduced by a factor of $10^{6}$. All the other
iterations are terminated when the residual norm is reduced by a factor of
$10^{4}$. The other assumptions are as the first part of the numerical
experiments.

Numerical results have been presented in Tables 2, 3 and 4. As seen, for large
values of $\nu$ and $\omega$ the ASSS preconditioner outperforms the other
preconditioners. However, for small values of $\nu$ and $\omega$ the ASSS
preconditioner is less effective than the others.

Finally, in Figure 1 the eigenvalue distribution of the matrices $\mathcal{B}$
and $\mathcal{P}_{\alpha}^{-1}\mathcal{B}$ with $\alpha^{*}$ (reported in
Table 1) have been displayed for $h=2^{-4}$ and different values of $\nu$ and
$\omega$. As we observe, the eigenvalues of the preconditioned matrix are
well-clustered in the circle with radius 1 and centered at $(1,0)$.

Table 2: Number of iterations of the methods along with the elapsed CPU time
(in parenthesis) for $h=2^{-5}$ and different values of $\nu$ and $\omega$.

| $\nu\setminus\omega$ | $10^{-4}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$ | $1$ | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$
---|---|---|---|---|---|---|---|---|---|---
IASSS | $10^{-2}$ | 54(0.12) | 54(0.11) | 54(0.10) | 54(0.10) | 54(0.10) | 53(0.09) | 45(0.07) | 40(0.06) | 51(0.06)
| $10^{-4}$ | 45(0.10) | 45(0.08) | 45(0.07) | 45(0.07) | 45(0.07) | 45(0.07) | 43(0.07) | 40(0.06) | 51(0.07)
| $10^{-6}$ | 40(0.09) | 40(0.06) | 40(0.06) | 40(0.06) | 40(0.06) | 40(0.06) | 40(0.06) | 42(0.06) | 51(0.07)
| $10^{-8}$ | 51(0.10) | 51(0.06) | 51(0.07) | 51(0.06) | 51(0.07) | 51(0.06) | 51(0.07) | 51(0.07) | 52(0.07)
IBAS | $10^{-2}$ | 38(0.11) | 38(0.07) | 38(0.10) | 38(0.06) | 38(0.06) | 24(0.04) | 473(0.60) | ${\dagger}$ | ${\dagger}$
| $10^{-4}$ | 33(0.09) | 33(0.05) | 33(0.05) | 33(0.05) | 33(0.06) | 33(0.05) | 39(0.06) | ${\dagger}$ | ${\dagger}$
| $10^{-6}$ | 33(0.08) | 33(0.05) | 33(0.05) | 33(0.06) | 33(0.05) | 33(0.05) | 33(0.05) | 57(0.07) | ${\dagger}$
| $10^{-8}$ | 38(0.07) | 38(0.06) | 38(0.06) | 38(0.06) | 38(0.06) | 38(0.06) | 38(0.06) | 38(0.06) | 68(0.08)
P-ASSS | $10^{-2}$ | 29(0.13) | 29(0.09) | 29(0.09) | 29(0.09) | 29(0.09) | 29(0.08) | 28(0.07) | 21(0.05) | 20(0.04)
| $10^{-4}$ | 28(0.11) | 28(0.08) | 28(0.07) | 28(0.07) | 28(0.08) | 28(0.07) | 28(0.07) | 21(0.04) | 20(0.04)
| $10^{-6}$ | 21(0.09) | 21(0.05) | 21(0.05) | 21(0.05) | 21(0.05) | 21(0.05) | 21(0.05) | 21(0.05) | 20(0.04)
| $10^{-8}$ | 20(0.09) | 20(0.04) | 20(0.04) | 20(0.04) | 20(0.04) | 20(0.04) | 20(0.04) | 20(0.04) | 19(0.04)
P-BAS | $10^{-2}$ | 18(0.12) | 19(0.08) | 20(0.05) | 20(0.05) | 20(0.05) | 17(0.06) | 26(0.07) | 48(0.12) | 36(0.06)
| $10^{-4}$ | 19(0.12) | 20(0.07) | 21(0.05) | 21(0.05) | 22(0.08) | 21(0.04) | 19(0.04) | 47(0.10) | 40(0.07)
| $10^{-6}$ | 18(0.11) | 19(0.06) | 20(0.04) | 21(0.06) | 21(0.08) | 22(0.04) | 22(0.03) | 28(0.06) | 42(0.08)
| $10^{-8}$ | 17(0.13) | 18(0.04) | 19(0.03) | 20(0.04) | 21(0.04) | 21(0.04) | 21(0.03) | 21(0.04) | 26(0.04)
P-PRESB | $10^{-2}$ | 7(0.20) | 7(0.05) | 7(0.11) | 8(0.10) | 10(0.16) | 23(0.75) | 106(8.80) | ${\ddagger}$ | ${\ddagger}$
| $10^{-4}$ | 8(0.18) | 8(0.04) | 8(0.07) | 8(0.09) | 9(0.09) | 12(0.20) | 60(3.06) | ${\ddagger}$ | ${\ddagger}$
| $10^{-6}$ | 8(0.20) | 8(0.05) | 8(0.04) | 8(0.07) | 8(0.06) | 9(0.09) | 12(0.18) | 89(9.78) | ${\ddagger}$
| $10^{-8}$ | 8(0.19) | 8(0.04) | 8(0.03) | 8(0.04) | 8(0.07) | 8(0.06) | 8(0.05) | 11(0.17) | 60(3.03)
P-BD | $10^{-2}$ | 14(0.15) | 14(0.07) | 14(0.07) | 16(0.07) | 22(0.08) | 54(0.20) | 176(1.27) | 132(0.72) | 26(0.06)
| $10^{-4}$ | 16(0.14) | 16(0.07) | 18(0.06) | 22(0.09) | 42(0.15) | 188(1.42) | 481(25.68) | 351(5.25) | 36(0.07)
| $10^{-6}$ | 15(0.15) | 15(0.06) | 15(0.04) | 16(0.04) | 30(0.10) | 138(0.77) | 345(5.21) | 298(3.49) | 34(0.07)
| $10^{-8}$ | 15(0.13) | 15(0.06) | 15(0.04) | 14(0.03) | 11(0.03) | 22(0.06) | 34(0.08) | 34(0.07) | 20(0.04)

Table 3: Number of iterations of the methods along with the elapsed CPU time
(in parenthesis) for $h=2^{-6}$ and different values of $\nu$ and $\omega$.

| $\nu\setminus\omega$ | $10^{-4}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$ | $1$ | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$
---|---|---|---|---|---|---|---|---|---|---
IASSS | $10^{-2}$ | 56(0.59) | 56(0.55) | 56(0.55) | 56(0.57) | 56(0.57) | 55(0.54) | 50(0.42) | 40(0.31) | 51(0.29)
| $10^{-4}$ | 50(0.45) | 50(0.42) | 50(0.41) | 50(0.41) | 50(0.41) | 50(0.41) | 48(0.39) | 40(0.32) | 51(0.30)
| $10^{-6}$ | 40(0.34) | 40(0.34) | 40(0.31) | 40(0.31) | 40(0.30) | 40(0.31) | 40(0.31) | 42(0.30) | 51(0.29)
| $10^{-8}$ | 51(0.32) | 51(0.29) | 51(0.30) | 51(0.29) | 51(0.30) | 51(0.30) | 51(0.29) | 51(0.31) | 52(0.30)
IBAS | $10^{-2}$ | 38(0.36) | 38(0.34) | 38(0.33) | 38(0.32) | 38(0.33) | 24(0.23) | 467(2.52) | ${\dagger}$ | ${\dagger}$
| $10^{-4}$ | 35(0.27) | 35(0.25) | 35(0.25) | 33(0.25) | 35(0.25) | 35(0.25) | 39(0.28) | ${\dagger}$ | ${\dagger}$
| $10^{-6}$ | 33(0.23) | 33(0.21) | 33(0.20) | 33(0.20) | 33(0.21) | 33(0.20) | 33(0.20) | 56(0.31) | ${\dagger}$
| $10^{-8}$ | 38(0.23) | 38(0.21) | 38(0.22) | 38(0.21) | 38(0.21) | 38(0.20) | 38(0.20) | 38(0.20) | 65(0.28)
P-ASSS | $10^{-2}$ | 28(0.44) | 28(0.41) | 28(0.40) | 28(0.38) | 28(0.38) | 29(0.42) | 29(0.35) | 21(0.25) | 21(0.16)
| $10^{-4}$ | 29(0.41) | 29(0.39) | 29(0.38) | 29(0.39) | 29(0.38) | 29(0.38) | 29(0.36) | 25(0.24) | 21(0.16)
| $10^{-6}$ | 25(0.30) | 25(0.26) | 25(0.27) | 25(0.26) | 25(0.25) | 25(0.24) | 25(0.25) | 25(0.23) | 21(0.16)
| $10^{-8}$ | 21(0.25) | 21(0.18) | 21(0.17) | 21(0.17) | 21(0.20) | 21(0.17) | 21(0.18) | 21(0.16) | 21(0.16)
P-BAS | $10^{-2}$ | 18(0.27) | 19(0.31) | 20(0.32) | 20(0.32) | 20(0.32) | 17(0.25) | 26(0.30) | 49(0.45) | 44(0.35)
| $10^{-4}$ | 20(0.26) | 20(0.25) | 21(0.26) | 22(0.29) | 22(0.25) | 22(0.27) | 20(0.20) | 47(0.44) | 48(0.35)
| $10^{-6}$ | 18(0.17) | 19(0.19) | 20(0.19) | 21(0.23) | 21(0.21) | 21(0.19) | 22(0.21) | 28(0.21) | 47(0.37)
| $10^{-8}$ | 17(0.13) | 19(0.16) | 20(0.16) | 20(0.17) | 21(0.17) | 21(0.19) | 22(0.18) | 22(0.20) | 28(0.20)
P-PRESB | $10^{-2}$ | 7(0.31) | 7(0.20) | 7(0.29) | 8(0.40) | 10(0.62) | 23(2.92) | 107(33.19) | ${\ddagger}$ | ${\ddagger}$
| $10^{-4}$ | 8(0.28) | 8(0.17) | 8(0.20) | 8(0.26) | 9(0.32) | 12(0.76) | 60(10.00) | ${\ddagger}$ | ${\ddagger}$
| $10^{-6}$ | 8(0.24) | 8(0.14) | 8(0.12) | 8(0.19) | 8(0.21) | 9(0.24) | 12(0.55) | 91(14.59) | ${\ddagger}$
| $10^{-8}$ | 8(0.22) | 8(0.10) | 8(0.11) | 8(0.11) | 8(0.17) | 8(0.18) | 9(0.17) | 12(0.14) | 81(11.82)
P-BD | $10^{-2}$ | 14(0.33) | 14(0.25) | 14(0.21) | 16(0.23) | 22(0.36) | 54(0.81) | 182(4.42) | 174(3.86) | 84(1.16)
| $10^{-4}$ | 16(0.30) | 16(0.21) | 16(0.20) | 22(0.26) | 42(0.54) | 195(5.09) | ${\ddagger}$ | ${\ddagger}$ | 154(3.02)
| $10^{-6}$ | 16(0.23) | 16(0.15) | 16(0.14) | 17(0.16) | 38(0.35) | 185(4.25) | ${\ddagger}$ | ${\ddagger}$ | 168(3.31)
| $10^{-8}$ | 15(0.18) | 15(0.12) | 15(0.10) | 15(0.09) | 22(0.17) | 85(1.07) | 153(2.97) | 166(3.27) | 116(1.80)

Table 4: Number of iterations of the methods along with the elapsed CPU time
(in parenthesis) for $h=2^{-7}$ and different values of $\nu$ and $\omega$.

| $\nu\setminus\omega$ | $10^{-4}$ | $10^{-3}$ | $10^{-2}$ | $10^{-1}$ | $1$ | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$
---|---|---|---|---|---|---|---|---|---|---
IASSS | $10^{-2}$ | 57(3.43) | 57(3.30) | 57(3.26) | 57(3.35) | 57(3.31) | 56(2.95) | 52(2.45) | 44(1.70) | 51(1.41)
| $10^{-4}$ | 52(2.46) | 52(2.36) | 52(2.43) | 52(2.39) | 52(2.43) | 52(2.56) | 52(2.48) | 44(1.78) | 51(1.42)
| $10^{-6}$ | 44(1.78) | 44(1.65) | 44(1.67) | 44(1.78) | 44(1.78) | 44(1.79) | 44(1.73) | 43(1.60) | 51(1.44)
| $10^{-8}$ | 51(1.53) | 51(1.49) | 51(1.51) | 51(1.42) | 51(1.43) | 51(1.45) | 51(1.50) | 51(1.45) | 52(1.35)
IBAS | $10^{-2}$ | ${\dagger}$ | ${\dagger}$ | ${\dagger}$ | ${\dagger}$ | 48(2.30) | 25(1.32) | 465(14.04) | ${\dagger}$ | ${\dagger}$
| $10^{-4}$ | 36(1.43) | 36(1.43) | 36(1.43) | 36(1.46) | 36(1.40) | 36(1.42) | 39(1.70) | ${\dagger}$ | ${\dagger}$
| $10^{-6}$ | 33(1.09) | 33(1.10) | 33(1.08) | 33(1.10) | 33(1.10) | 33(1.14) | 33(1.13) | 56(1.78) | ${\dagger}$
| $10^{-8}$ | 38(0.98) | 38(0.98) | 38(0.97) | 38(0.97) | 38(0.96) | 38(0.96) | 38(0.97) | 38(0.97) | 64(1.48)
P-ASSS | $10^{-2}$ | 28(2.40) | 28(2.40) | 28(2.40) | 28(2.40) | 28(2.40) | 29(2.42) | 30(2.42) | 28(1.40) | 25(0.92)
| $10^{-4}$ | 30(2.38) | 30(2.42) | 30(2.44) | 30(2.38) | 30(2.54) | 30(2.37) | 29(2.23) | 28(1.51) | 25(0.98)
| $10^{-6}$ | 28(1.55) | 28(1.59) | 28(1.57) | 28(1.52) | 28(1.50) | 28(1.49) | 28(1.48) | 27(1.26) | 25(0.95)
| $10^{-8}$ | 25(0.95) | 25(0.94) | 25(0.94) | 25(0.94) | 25(0.94) | 25(0.94) | 25(0.92) | 25(0.90) | 24(0.79)
P-BAS | $10^{-2}$ | 18(1.77) | 19(1.84) | 20(1.87) | 20(1.88) | 20(1.83) | 17(1.61) | 26(1.77) | 49(2.33) | 46(1.47)
| $10^{-4}$ | 20(1.35) | 21(1.42) | 22(1.41) | 22(1.50) | 22(1.47) | 22(1.54) | 20(1.30) | 47(2.15) | 50(1.71)
| $10^{-6}$ | 18(0.76) | 19(0.77) | 20(0.80) | 21(0.90) | 21(0.81) | 21(0.87) | 22(0.96) | 28(1.16) | 49(1.67)
| $10^{-8}$ | 18(0.48) | 19(0.52) | 20(0.56) | 30(0.53) | 21(0.59) | 22(0.62) | 22(0.64) | 22(0.59) | 28(0.81)
P-PRESB | $10^{-2}$ | 7(1.15) | 7(1.09) | 7(1.92) | 8(2.48) | 10(3.79) | 23(18.42) | ${\ddagger}$ | ${\ddagger}$ | ${\ddagger}$
| $10^{-4}$ | 8(1.04) | 8(0.92) | 8(1.25) | 8(1.62) | 9(1.73) | 12(4.16) | 60(58.46) | ${\ddagger}$ | ${\ddagger}$
| $10^{-6}$ | 8(0.68) | 8(0.60) | 8(0.58) | 8(0.95) | 8(0.97) | 9(1.34) | 12(2.56) | 92(70.74) | ${\ddagger}$
| $10^{-8}$ | 8(0.51) | 8(0.41) | 8(0.40) | 8(0.40) | 8(0.65) | 8(0.62) | 9(0.81) | 12(1.66) | 87(47.56)
P-BD | $10^{-2}$ | 14(1.45) | 14(1.33) | 14(1.35) | 16(1.50) | 22(1.90) | 54(4.12) | 182(17.19) | 188(15.20) | 142(9.10)
| $10^{-4}$ | 16(1.21) | 16(1.24) | 18(1.27) | 22(1.51) | 42(2.76) | 199(19.97) | ${\ddagger}$ | ${\ddagger}$ | ${\ddagger}$
| $10^{-6}$ | 16(0.75) | 16(0.67) | 16(0.66) | 20(0.79) | 40(1.64) | 203(17.87) | ${\ddagger}$ | ${\ddagger}$ | ${\ddagger}$
| $10^{-8}$ | 15(0.47) | 15(0.41) | 15(0.40) | 16(0.42) | 30(0.89) | 151(10.12) | ${\ddagger}$ | ${\ddagger}$ | ${\dagger}$

Figure 1: Eigenvalue distribution of the matrices $\mathcal{B}$ and
$\mathcal{P}_{\alpha}^{-1}\mathcal{B}$ for $h=2^{-4}$ with different values of
$\nu$ and $\omega$. The value of $\alpha$ is
$\alpha^{*}=\sqrt{\mu_{\min}\mu_{\max}},$ where $\mu_{\min}$ and $\mu_{\max}$
are the smallest and largest eigenvalues of the matrix $M$, respectively.

##  6 Conclusion

We have presented a method called alternating symmetric positive definite and
scaled symmetric positive semidefinite splitting (ASSS) method, for solving
the system arisen from finite element discretization of a distributed optimal
control problem with time-periodic parabolic equations. We have proved that
the method is unconditionally convergent. We have compared the numerical
results of the ASSS method and the corresponding induced preconditioner with
those of the BAS iteration method. We also compared the numerical results of
the ASSS preconditioner to those of the PRESB preconditioner. Numerical
results showed that the proposed  method has some advantages over the two
other tested methods.

## Acknowledgments

The author would like to thank the referees for their careful reading of the
paper and giving several valuable comments and suggestions.

## References

  * [1] O. Axelsson, D. Lukas, Preconditioning methods for eddy-current optimally controlled time-harmonic electromagnetic problems, J. Numer. Math. 27 (2019) 1-21.
  * [2] O. Axelsson, D.K. Salkuyeh, A new version of a preconditioning method for certain two-by-two block matrices with square blocks, BIT Numer. Math. 59 (2018) 321–-342.
  * [3] O. Axelsson, M. Neytcheva, B. Ahmad, A comparison of iterative methods to solve complex valued linear algebraic systems, Numer. Algor. 66 (2014) 811–841.
  * [4] O. Axelsson, M. Neytcheva and Z.-Z. Liang, Parallel solution methods and preconditioners for evolution equations, Math. Model. Anal. 23 (2018) 287–308.
  * [5] O. Axelsson, P. Boyanova, M. Kronbichler, M. Neytcheva, X. Wu, Numerical and computational efficiency of solvers for two-phase problems, Comput. Math. Appl. 65 (2013) 301–314.
  * [6] Z.-Z. Bai, M. Benzi and F. Chen, Modified HSS iteration methods for a class of complex symmetric linear systems, Computing 87 (2010) 93–111.
  * [7] M. Benzi and G.H. Golub, A preconditioner for generalized saddle point problems, SIAM J. Matrix Anal. Appl. 26 (2004) 20–41.
  * [8] Y. Cao, A block positive-semidefinite splitting preconditioner for generalized saddle point linear systems, J. Comput. Appl. Math. 374 (2020) 112787.
  * [9] F. Chen, On choices of iteration parameter in HSS method, Appl. Math. Comput. 271 (2015) 832–837.
  * [10] H. Elman, D. Silvester, A.J. Wathen, Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics, Oxford University Press, 2005.
  * [11] Y.-M. Huang, A practical formula for computing optimal parameters in the HSS iteration methods, J. Comput. Appl. Math. 255 (2014) 142–149.
  * [12] M. Kollmann, M. Kolmbauer, A preconditioned MinRes solver for time-periodic parabolic optimal control problems, Numer. Linear Algebra Appl. 20 (2012) 761–784.
  * [13] W. Krendl, V. Simoncini and W. Zulehner, Stability estimates and structural spectral properties of saddle point problems, Numer. Math. 124 (2013) 183–213.
  * [14] J. Lions, Optimal Control of Systems, Springer, New York, 1968.
  * [15] H. Mirchi, D.K. Salkuyeh, A new iterative method for solving the systems arisen from finite element discretization of a time-harmonic parabolic optimal control problems, Math. Comput. Simul. 185 (2021) 771–782.
  * [16] T. Rees, H.S. Dollar and A.J. Wathen, Optimal solvers for PDE-constrained optimization, SIAM J. Sci. Comput. 32 (2010) 271–298.
  * [17] Z.-R. Ren, Y. Cao, An alternating positive-semidefinite splitting preconditioner for saddle point problems from time-harmonic eddy current models, IMA J. Numer. Anal 36 (2016) 922–946.
  * [18] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM J. Sci. Comput. 14 (1993) 461–469.
  * [19] Y. Saad, Iterative methods for sparse linear systems, Second Edition, SIAM, 2003.
  * [20] Y. Saad, M.H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869.
  * [21] A.J. Wathen, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA J. Numer. Anal. 7 (1987) 449–457.
  * [22] A.-L. Yang, Scaled norm minimization method for computing the parameters of the HSS and the two-parameter HSS preconditioners, Numer. Linear Algebra Appl. 25 (2018) e2169.
  * [23] M.-L. Zeng, Respectively scaled splitting iteration method for a class of block 4-by-4 linear systems from eddy current electromagnetic problems, Japan J. Indust. Appl. Math. (2020). https://doi.org/10.1007/s13160-020-00446-8.
  * [24] Z. Zheng, G.-F. Zhang and M.-Z. Zhu, A block alternating splitting iteration method for a class of block two-by-two complex linear systems, J. Comput. Math. Appl. 288 (2015) 203–214.
  * [25] Z. Zheng, G.-F. Zhang., M.-Z. Zhu, A note on preconditioners for complex linear systems arising from PDE-constrained optimization problems, Appl. Math. Lett. 61 (2016) 114–121.