---
abstract: 'We consider the solution of complex symmetric shifted linear systems. Such systems arise in large-scale electronic structure simulations and there is a strong need for the fast solution of the systems. With the aim of solving the systems efficiently, we consider a special case of the QMR method for non-Hermitian shifted linear systems and propose its weighted quasi-minimal residual approach. A numerical algorithm, referred to as shifted QMR\_SYM($B$), is given by the choice of a particularly cost-effective weight. Numerical examples are presented to show the performance of the shifted QMR\_SYM($B$) method.'
author:
- 'T. Sogabe[^1]'
- 'T. Hoshi[^2]'
- 'S.-L. Zhang$^\dagger$'
- 'T. Fujiwara[^3]'
title: 'On a weighted quasi-residual minimization strategy of the QMR method for solving complex symmetric shifted linear systems[^4]'
---

Introduction
============

In this paper we consider the solution of complex symmetric shifted linear systems of the form $$\begin{aligned}
  (A+\sigma_\ell I){\mbox{\boldmath $x$}}^{(\ell)}={\mbox{\boldmath $b$}}, \quad   \ell=1,2,\ldots,m,  \label{CS_SLS}\end{aligned}$$ where $A({\sigma_\ell}):=A+\sigma_\ell I$ are nonsingular $N$-by-$N$ complex symmetric sparse matrices, i.e., $A(\sigma_\ell)=A(\sigma_\ell)^{\textrm{\scriptsize{T}}}\neq \overline{A}(\sigma_\ell)^{\textrm{\scriptsize{T}}}$, with scalar shifts $\sigma_\ell \in \mathbf{C}$, $I$ is the $N$-by-$N$ identity matrix, and ${\mbox{\boldmath $x$}}^{(\ell)}, {\mbox{\boldmath $b$}}$ are complex vectors of length $N$. The above systems arise in large-scale electronic structure simulations [@Shifted_COCG], and there is a strong need for the fast solution of the systems.

Since the given shifted linear systems (\[CS\_SLS\]) are a set of sparse linear systems, it is natural to use Krylov subspace methods, and moreover since the coefficient matrices are complex symmetric, one of the simplest ways to solve the shifted linear systems is applying (preconditioned) Krylov subspace methods for solving complex symmetric linear systems such as the COCG method [@COCG], the COCR method [@COCR], and the QMR\_SYM method [@QMR_SYM] to all of the shifted linear systems (\[CS\_SLS\]). On the other hand, denoting the $n$-dimensional Krylov subspace with respect to $A$ and ${\mbox{\boldmath $b$}}$ as $K_n(A,{\mbox{\boldmath $b$}}):={\textrm}{span}\{{\mbox{\boldmath $b$}},A{\mbox{\boldmath $b$}},\ldots,A^{n-1}{\mbox{\boldmath $b$}}\}$, we observe that $$\begin{aligned}
  K_n(A,{\mbox{\boldmath $b$}})=K_n(A({\sigma_\ell}),{\mbox{\boldmath $b$}}).       \label{KS}\end{aligned}$$ This implies that once basis vectors are generated for one of the Krylov subspaces $K_n(A({\sigma_\ell}),{\mbox{\boldmath $b$}})$, these basis vectors can be used to solve all the shifted linear systems. In other words, there is no need to generate all Krylov subspaces $K_n(A(\sigma_\ell),{\mbox{\boldmath $b$}})$, and thus computational costs involving the basis generation, e.g., matrix-vector multiplications, are saved. Here we give a concrete example: if we apply the conjugate orthogonal conjugate gradient (COCG) method to all the linear systems (\[CS\_SLS\]), then bases for $K_n(A(\sigma_\ell),{\mbox{\boldmath $b$}})$ are generated for $\ell=1,2,\ldots,m$. On the other hand, if we apply the COCG method to just one of the shifted linear systems (\[CS\_SLS\]) (referred to as the “$seed\ system$”), then the Krylov basis vectors are generated from the seed system and these vectors are used to solve the rest of the shifted linear systems.

Based on the observation (\[KS\]), the shifted COCG method [@Shifted_COCG] has been recently proposed for solving complex symmetric shifted linear systems. The feature of the shifted COCG method is that the method performs COCG on a seed system and makes it possible to complete COCG for all shifted linear systems without further matrix-vector multiplications. The feature is completely different from some of the other well-known shifted linear solvers such as the shifted BiCGStab($\ell$) method [@Shifted_BiCGSSTABL] or the restarted shifted GMRES method [@Shifted_GMRES] since these perform BiCGStab($\ell$) (or GMRES) on a seed system but a different method on the shifted linear systems in order to keep the residuals colinear. The feature of the shifted COCG method plays a very important role in the seed switching technique [@Seed_switch] that avoids a minor problem of the shifted COCG method: one requires the choice of a seed system and an unsuitable choice may lead to the undesirable result that some shifted linear systems remain unsolved.

There is another approach to solving the shifted linear systems (\[CS\_SLS\]). That is the use of Krylov subspace methods for non-Hermitian shifted linear systems such as the shifted BiCGStab($\ell$) method [@Shifted_BiCGSSTABL], the shifted (TF)QMR method [@Freund], the restarted shifted FOM method [@Shifted_FOM], and the restarted shifted GMRES method [@Shifted_GMRES], see also, e.g., [@Simoncini]. We readily see that the relation (\[KS\]) holds not only for complex symmetric matrices but also for non-Hermitian matrices, and these methods are based on the use of this shift-invariant relation. Therefore, this can be a good approach. However, since these methods do not exploit the property of complex symmetric matrices, their computational costs can be more expensive than that of the shifted COCG method.

In this paper we consider the shifted QMR\_SYM method that is a special case of the QMR method for non-Hermitian shifted linear systems [@Freund] and clarify the most time consuming part of it for a large number of shifted linear systems. Then, in order to reduce the cost, we propose a weighted quasi-minimal residual (WQMR) approach and propose a specific weight. We experimentally show the practical efficiency of the resulting algorithm, referred to as shifted QMR\_SYM($B$), when the number of shifted linear systems is large enough.

The present paper is organized as follows: in the next section, shifted QMR\_SYM is described in order to specify the most time consuming part for a large number of shifted linear systems. In Section 3, we propose a WQMR approach with a specific weight for reducing the cost of the most time consuming part. The resulting algorithm, shifted QMR\_SYM($B$), and its properties are given. In Section 4, some results of numerical examples from electronic structure simulations are shown to see the performance of the shifted QMR\_SYM($B$) method. Finally, some concluding remarks are made in Section 5.

Throughout this paper, unless otherwise stated, all vectors and matrices are assumed to be complex. $\overline{M}$, $M^{\textrm{\scriptsize{T}}}$, $M^{\textrm{\scriptsize{H}}}=\overline{M}^{\textrm{\scriptsize{T}}}$ denote the complex conjugate, transpose, and Hermitian transpose matrix of the matrix $M$. $\|{\mbox{\boldmath $v$}}\|_W$ denotes the $W$-norm written as $({\mbox{\boldmath $v$}}^{\textrm{\scriptsize{H}}}W {\mbox{\boldmath $v$}})^{1/2}$, where $W$ is Hermitian positive definite.

The QMR\_SYM method for solving complex symmetric shifted linear systems
========================================================================

In this section, the shifted QMR\_SYM method and its properties for solving complex symmetric shifted linear systems are introduced.

The QMR method for shifted linear systems was first formulated in [@Freund] for the case of a general non-Hermitian matrix. Therefore, by simplifying the non-Hermitian Lanczos process [@Lanczos], as is known from other papers such as [@QMR_SYM; @COCG], a shifted simplified QMR method, shifted QMR\_SYM, is readily obtained for the case of a complex symmetric matrix. Its algorithm is given next.\
\
$$\begin{aligned}
&&{\hspace{-7mm}}\\[-15mm]
&&{\hspace{-7mm}}{\mbox{\boldmath $x$}}_{0}^{(\ell)}={\mbox{\boldmath $p$}}_{-1}^{(\ell)}=
      {\mbox{\boldmath $p$}}_{0}^{(\ell)}=\mathbf{0}, \ 
      {\mbox{\boldmath $v$}}_{1}={\mbox{\boldmath $b$}}/({\mbox{\boldmath $b$}}^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}})^{1/2},\ g_1^{(\ell)}=({\mbox{\boldmath $b$}}^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}})^{1/2},\\
&&{\hspace{-7mm}}{\textrm}{{\bf for}}\ n=1,2, \ldots\ {\textrm}{{\bf do:}}\\
&&{\hspace{-7mm}}\quad {\textrm}{(The complex symmetric Lanczos process)}\nonumber\\
&&{\hspace{-7mm}}\quad \alpha_n ={\mbox{\boldmath $v$}}_n^{\textrm{\scriptsize{T}}}A{\mbox{\boldmath $v$}}_n,                           \nonumber \\
&&{\hspace{-7mm}}\quad \widetilde{{\mbox{\boldmath $v$}}}_{n+1}=
        A{\mbox{\boldmath $v$}}_n-\alpha_n{\mbox{\boldmath $v$}}_n-\beta_{n-1}{\mbox{\boldmath $v$}}_{n-1},              \nonumber \\
&&{\hspace{-7mm}}\quad \beta_{n} =(\widetilde{{\mbox{\boldmath $v$}}}_{n+1}^{\textrm{\scriptsize{T}}}\widetilde{{\mbox{\boldmath $v$}}}_{n+1})^{1/2},                 \nonumber \\
&&{\hspace{-7mm}}\quad {\mbox{\boldmath $v$}}_{n+1} = \widetilde{{\mbox{\boldmath $v$}}}_{n+1}/\beta_{n},          \nonumber \\
&&{\hspace{-7mm}}\quad t_{n-1,n}^{(\ell)}=\beta_{n-1},\ t_{n,n}^{(\ell)} =\alpha_{n}+\sigma_{\ell},\  t_{n+1,n}^{(\ell)}=\beta_{n},         \nonumber \\
&&{\hspace{-7mm}}\quad {\textrm}{(Solve least squares problems by Givens rotations)}                            \nonumber \\
&&{\hspace{-7mm}}\quad  {\textrm}{{\bf for}}\ \ell=1,2, \ldots, m\ {\textrm}{{\bf do:}}\\
&&{\hspace{-7mm}}\qquad {\textrm}{{\bf if}}\ \|{\mbox{\boldmath $r$}}_{n}^{(\ell)}\|_2/\|{\mbox{\boldmath $b$}}\|_2> \epsilon,\
            {\textrm}{{\bf then}}\\
&&{\hspace{-7mm}}\qquad \quad {\textrm}{{\bf for}}\ i=\max\{1,n-2\}, \ldots, n-1\ 
                   {\textrm}{{\bf do:}}\\
&&{\hspace{-7mm}}\qquad \qquad \pmatrix{
                t_{i,n}^{(\ell)} \cr
               t_{i+1,n}^{(\ell)}}
=
\pmatrix{
           c_{i}^{(\ell)}     & s_{i}^{(\ell)}  \cr
           -\overline{s}_{i}^{(\ell)} & c_{i}^{(\ell)} }
\pmatrix{
               t_{i,n}^{(\ell)} \cr
               t_{i+1,n}^{(\ell)} },      \nonumber\\
&&{\hspace{-7mm}}\qquad \quad {\textrm}{{\bf end}}             \nonumber\\
&&{\hspace{-7mm}}\qquad \quad c_n^{(\ell)}=\frac{|t_{n,n}^{(\ell)}|}{\sqrt{|t_{n,n}^{(\ell)}|^2+|t_{n+1,n}^{(\ell)}|^2}}, \nonumber \\
&&{\hspace{-7mm}}\qquad  \quad \overline{s}_n^{(\ell)}=\frac{t_{n+1,n}^{(\ell)}}{t_{n,n}^{(\ell)}}c_n^{(\ell)},   \nonumber \\ 
&&{\hspace{-7mm}}\qquad  \quad  t_{n,n}^{(\ell)}  =c_n^{(\ell)}t_{n,n}^{(\ell)}+s_n^{(\ell)}t_{n+1,n}^{(\ell)},  \nonumber \\
&&{\hspace{-7mm}}\qquad  \quad t_{n+1,n}^{(\ell)}=0,                        \nonumber \\
&&{\hspace{-7mm}}\qquad \quad \pmatrix{
                g_{n}^{(\ell)} \cr
                g_{n+1}^{(\ell)} }
 =
 \pmatrix{
           c_{n}^{(\ell)}     & s_{n}^{(\ell)}  \cr
           -\overline{s}_{n}^{(\ell)} & c_{n}^{(\ell)} }
\pmatrix{
               g_{n}^{(\ell)} \cr
               0 },   \nonumber\\   
&&{\hspace{-7mm}}\qquad \quad  {\textrm}{(Update approximate solutions ${\mbox{\boldmath $x$}}_n^{(\ell)}$)}\nonumber\\
&&{\hspace{-7mm}}\qquad \quad {\mbox{\boldmath $p$}}_{n}^{(\ell)}= {\mbox{\boldmath $v$}}_{n}-(t_{n-2,n}^{(\ell)}/t_{n-2,n-2}^{(\ell)}){\mbox{\boldmath $p$}}_{n-2}^{(\ell)}- (t_{n-1,n}^{(\ell)}/t_{n-1,n-1}^{(\ell)}){\mbox{\boldmath $p$}}_{n-1}^{(\ell)}, \nonumber\\
&&{\hspace{-7mm}}\qquad \quad {\mbox{\boldmath $x$}}_{n}^{(\ell)}={\mbox{\boldmath $x$}}_{n-1}^{(\ell)}+(g_n^{(\ell)}/t_{n,n}^{(\ell)}){\mbox{\boldmath $p$}}_n^{(\ell)},  \nonumber \\
&&{\hspace{-7mm}}\quad \quad {\textrm}{{\bf end if}}       \nonumber    \\
&&{\hspace{-7mm}}\quad {\textrm}{{\bf end}}                \nonumber    \\&&{\hspace{-7mm}}\qquad {\textrm}{{\bf if}}\ \|{\mbox{\boldmath $r$}}_{n}^{(\ell)}\|_2/\|{\mbox{\boldmath $b$}}\|_2\le \epsilon \ 
             {\textrm}{for all}\ \ell,\ {\textrm}{{\bf then} exit}.\\
&&{\hspace{-7mm}}{\textrm}{{\bf end}}                      \nonumber             \end{aligned}$$ Algorithm 2.1 can be regarded as a natural combination of the results given in [@QMR_SYM; @Freund].

In order to know that the numerical solution is accurate enough, one may need to compute the residual 2-norms. In that case, the following computation may be useful.

The 2-norms of the $n$th residuals of the approximate solutions ${\mbox{\boldmath $x$}}_n^{(\ell)}$ of the shifted QMR\_SYM method are given by $$\begin{aligned}
 \|{\mbox{\boldmath $r$}}_n^{(\ell)}\|_2=|g_{n+1}^{(\ell)}|\cdot\|{\mbox{\boldmath $w$}}_{n+1}^{(\ell)}\|_2\quad {\textrm}{for} \quad   \ell=1,2,\ldots,m, \end{aligned}$$ where ${\mbox{\boldmath $w$}}_{n+1}^{(\ell)}=-s_n^{(\ell)}{\mbox{\boldmath $w$}}_n^{(\ell)}+c_n^{(\ell)}{\mbox{\boldmath $v$}}_{n+1}$ and ${\mbox{\boldmath $w$}}_1^{(\ell)}={\mbox{\boldmath $v$}}_1.$

Proposition 2.1 is a result known to hold for the QMR method [@FrNa91]. Therefore, it also holds for the above specialized variant. The rest of this section describes some special properties of the shifted QMR\_SYM method.

Let $A\in \mathbf{R}^{N\times N}$ be real symmetric, $\sigma_\ell\in\mathbf{C}$ be complex shifts, and ${\mbox{\boldmath $b$}}\in\mathbf{R}^N$. Then the shifted QMR\_SYM method (Algorithm 2.1) enjoys the following properties:

-   all matrix-vector multiplications can be done in real arithmetic;

-   an approximate solution at $n$th iteration step for each $\ell$ has minimal residual 2-norms, i.e., ${\mbox{\boldmath $x$}}_n^{(\ell)}$’s are generated such that $\min\|{\mbox{\boldmath $r$}}_n^{(\ell)}\|_2$ over ${\mbox{\boldmath $x$}}_n^{(\ell)}\in K_n(A,{\mbox{\boldmath $b$}})$;

-   $ \|{\mbox{\boldmath $r$}}_n^{(\ell)}\|_2=|g_{n+1}^{(\ell)}|\  {\textrm}{for} \  \ell=1,2,\ldots,m, \  n\ge 0.$

The above properties are known results since the properties have been proved for each individual shift. See [@Fr90] for details. The properties of Proposition 2.2 may be very useful for large-scale electronic structure simulations [@Shifted_COCG] and a projection approach for eigenvalue problems [@Sakurai] since there are complex symmetric shifted linear systems to be solved efficiently under the assumption of Proposition 2.2.

An iterative method for solving complex symmetric shifted linear systems
========================================================================

In this section we consider complex symmetric shifted linear systems with a large number of shifts. For such systems, say $m\gg 1$, the most time consuming part of Algorithm 2.1 can be generating approximate solutions since the cost for the recurrences is $6Nm+3m$ per iteration step. In order to reduce the cost, in this section we propose a weighted quasi-minimal residual approach with a specific weight. In the next subsection we discuss the details of the approach. In Subsection 3.2 we give a specific weight to achieve the reduction of the cost.

A weighted quasi-minimal residual approach
------------------------------------------

Before the formulation of a Weighted Quasi-Minimal Residual (WQMR) approach, let us recall the complex symmetric Lanczos process (see, e.g., Algorithm 2.1 in [@QMR_SYM]).\
\
$$\begin{aligned}
&&{\hspace{-7mm}}{\textrm}{set}\ \beta_0=0,\ {\mbox{\boldmath $v$}}_0=\mathbf{0},\ {\mbox{\boldmath $r$}}_0\ne \mathbf{0}\in \mathbf{C}^N,\ \nonumber \\
&&{\hspace{-7mm}}{\textrm}{set}\ {\mbox{\boldmath $v$}}_1={\mbox{\boldmath $r$}}_0/({\mbox{\boldmath $r$}}_0^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $r$}}_0)^{1/2},         \nonumber \\
&&{\hspace{-7mm}}{\textrm}{for $n=1,2,\ldots, m-1$ do:}              \nonumber \\
&&{\hspace{-7mm}}\quad \alpha_n ={\mbox{\boldmath $v$}}_n^{\textrm{\scriptsize{T}}}A{\mbox{\boldmath $v$}}_n,                      \nonumber \\
&&{\hspace{-7mm}}\quad \widetilde{{\mbox{\boldmath $v$}}}_{n+1}=
        A{\mbox{\boldmath $v$}}_n-\alpha_n{\mbox{\boldmath $v$}}_n-\beta_{n-1}{\mbox{\boldmath $v$}}_{n-1},              \nonumber \\
&&{\hspace{-7mm}}\quad \beta_{n} =(\widetilde{{\mbox{\boldmath $v$}}}_{n+1}^{\textrm{\scriptsize{T}}}\widetilde{{\mbox{\boldmath $v$}}}_{n+1})^{1/2},   \nonumber \\
&&{\hspace{-7mm}}\quad {\mbox{\boldmath $v$}}_{n+1} = \widetilde{{\mbox{\boldmath $v$}}}_{n+1}/\beta_{n}.               \nonumber \\
&&{\hspace{-7mm}}{\textrm}{end}                                   \nonumber\end{aligned}$$ The matrix form of Algorithm 3.1 is known as follows: let $T_{n+1,n}$ and $T_n$ be the $(n+1)\times n$ and $n\times n$ tridiagonal matrices whose entries are recurrence coefficients of the complex symmetric Lanczos process, which are given by $$\begin{aligned}
T_{n+1,n}:=
\pmatrix{
    \alpha_1& \beta_1 &        &                \cr
    \beta_1 & \alpha_2& \ddots &                \cr
            & \ddots  & \ddots & \beta_{n-1}    \cr
            &         & \beta_{n-1} & \alpha_n  \cr
            &         &        & \beta_{n}},    \quad 
T_{n}:=
\pmatrix{
    \alpha_1& \beta_1 &        &                \cr
    \beta_1 & \alpha_2& \ddots &                \cr
            & \ddots  & \ddots & \beta_{n-1}    \cr
            &         & \beta_{n-1} & \alpha_n },\end{aligned}$$ and let $V_n$ be the $N\times n$ matrix with the Lanczos vectors as columns, i.e., $V_n:=({\mbox{\boldmath $v$}}_1,{\mbox{\boldmath $v$}}_2,\ldots,{\mbox{\boldmath $v$}}_n)$. Then from Algorithm 3.1, it follows that $$\begin{aligned}
&&\hspace{4mm} AV_n = V_{n+1}T_{n+1,n}  
  = V_n T_n +\beta_{n}{\mbox{\boldmath $v$}}_{n+1}{\mbox{\boldmath $e$}}^{{\textrm{\scriptsize{T}}}}_n, \label{comp_sym_Lanczos_Mat}\end{aligned}$$ where ${\mbox{\boldmath $e$}}_n=(0,0,\ldots,1)^{\textrm{\scriptsize{T}}}\in \mathbf{R}^n$.

Now we are ready to describe the WQMR approach. Let ${\mbox{\boldmath $x$}}_n^{(\ell)}$ be approximate solutions at the $n$th iteration step for the systems (\[CS\_SLS\]), which are given by $${\mbox{\boldmath $x$}}_n^{(\ell)} = V_n{\mbox{\boldmath $y$}}_n^{(\ell)},\quad \ell=1,2,\ldots,m, \label{sol}$$ where ${\mbox{\boldmath $y$}}_n^{(\ell)} \in \mathbf{C}^n$. Then, from the definition of residual vectors ${\mbox{\boldmath $r$}}_n^{(\ell)}:={\mbox{\boldmath $b$}}-(A+\sigma_\ell I){\mbox{\boldmath $x$}}_n^{(\ell)}$, the update formulas (\[sol\]), and the matrix form of the complex symmetric Lanczos process (\[comp\_sym\_Lanczos\_Mat\]), we readily obtain $$\begin{aligned}
{\mbox{\boldmath $r$}}_n^{(\ell)}  = V_{n+1}\Bigl(g_1{\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $y$}}_n^{(\ell)}\Bigr),\ \ {\textrm}{where}\ \  T_{n+1,n}^{(\ell)}:= T_{n+1,n}+ \sigma_\ell
            \pmatrix{
              I_n \cr
              \mathbf{0}^{\textrm{\scriptsize{T}}}\cr 
            }.                            \label{QMR_residual}\end{aligned}$$ Here, ${\mbox{\boldmath $e$}}_1$ is the first unit vector written by ${\mbox{\boldmath $e$}}_1$=$(1,0,\ldots,0)^{\textrm{\scriptsize{T}}}$ and $g_1=({\mbox{\boldmath $b$}}^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}})^{1/2}$. It is natural to determine ${\mbox{\boldmath $y$}}_n^{(\ell)}$ such that all residual 2-norms $\|{\mbox{\boldmath $r$}}_n^{(\ell)}\|_2$ are minimized. However, this choice for ${\mbox{\boldmath $y$}}_n^{(\ell)}$ is impractical due to large computational costs. Hence, an alternative approach is given, i.e., the vectors ${\mbox{\boldmath $y$}}_n^{(\ell)}$ are determined by solving the following weighted least squares problems: $$\begin{aligned}
{\mbox{\boldmath $y$}}_n^{(\ell)}  = \arg\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\| g_1{\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\Bigr\|_{W_{n+1}^{\textrm{\scriptsize{H}}}W_{n+1}^{}},  
                                           \label{WLSP}\end{aligned}$$ where $W_{n+1}$ is an $(n+1)$-by-$(n+1)$ nonsingular matrix. Thus $W_{n+1}^{\textrm{\scriptsize{H}}}W_{n+1}$ can be used as a weight since it is Hermitian positive definite. One of the simplest choices for $W_{n+1}$ is the identity matrix. In this case, from $W_{n+1}=I_{n+1}$ we have $$\begin{aligned}
{\mbox{\boldmath $y$}}_n^{(\ell)}  = \arg\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\| g_1 {\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\Bigr\|_2.  \label{WLSP_2norm}\end{aligned}$$ The vector that is minimized is called quasi-residual. Algorithm 2.1 is obtained by solving (\[WLSP\_2norm\]) using Givens rotations, see, e.g., [@Mat_Comp p.215].

A slightly generalized choice proposed in [@FrNa91] is $W_{n+1}=\Omega_{n+1}:={\textrm}{diag}(\omega_1,\omega_2,\ldots,\\ \omega_{n+1})$ with $\omega_i>0$ for all $i$. Then, it follows that $$\begin{aligned}
{\mbox{\boldmath $y$}}_n^{(\ell)}  = \arg \min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\| \omega_1g_1{\mbox{\boldmath $e$}}_1- \Omega_{n+1} T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\Bigr\|_2.\end{aligned}$$ Of various possible choices for $\omega_i$, a natural one is $\omega_i=\|{\mbox{\boldmath $v$}}_i\|_2$ since $\Omega_{n+1}$ contains then the diagonal entries of the upper triangular matrix $R_{n+1}$ that is obtained by the $QR$ factorization of $V_{n+1}$. If we choose $W_{n+1}=R_{n+1}$, where $V_{n+1}=Q_{n+1}R_{n+1}$, then from (\[QMR\_residual\]) and (\[WLSP\]) we have $$\begin{aligned}
\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\|g_1{\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)} \Bigr\|_{R_{n+1}^{\textrm{\scriptsize{H}}}R_{n+1}^{}}
&=&\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\|g_1R_{n+1}{\mbox{\boldmath $e$}}_1- R_{n+1}T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)} \Bigr\|_2 \\
&=&\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\|Q_{n+1}R_{n+1}\bigl(g_1{\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\bigr)\Bigr\|_2 \\
&=&\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\|V_{n+1}\bigl(g_1{\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\bigr)\Bigr\|_2 \\
&=&\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \bigl\|{\mbox{\boldmath $r$}}_n^{(\ell)}\bigr\|_2. \\\end{aligned}$$ By solving the above weighted least squares problems, all residual 2-norms are minimized. Hence $W_{n+1}=\Omega_{n+1}$ is a rational choice. For each individual shift, the resulting algorithm is the same as that in [@QMR_SYM Algorithm 3.2].

A choice of the weight suitable for a large number of shifts
------------------------------------------------------------

In the previous subsection, we have described the WQMR approach and mentioned that the choice of the weight $W_{n+1}^{\textrm{\scriptsize{H}}}W_{n+1}$ with $W_{n+1}=I_{n+1}$ leads to the shifted QMR\_SYM method (Algorithm 2.1). Under the assumption of proposition 2.2, the shifted QMR\_SYM method is ideal in the sense of Faber-Manteuffel’s theorem [@Faber] since it enjoys minimal residual property and requires not long-term but short-term recurrences for updating approximate solutions, and thus one may think that there is no need to choose other possible weights. However, we will show in this subsection that even under the assumption of proposition 2.2 there is a practical weight for the WQMR approach. The motivation for the choice of the weight mainly comes from the freedom of the number $m$ of complex symmetric shifted linear systems.

Now we consider the computational costs of Algorithm 2.1 for a large number of complex symmetric shifted linear systems, i.e., $m\gg1$. For the case $m\gg1$ we readily see from Algorithm 2.1 that computing the recurrences for updating approximate solutions is the most time-consuming part due to a cost of $6Nm+3m$ per iteration step. Hence we will now consider a weight to reduce the computational cost for the recurrences of ${\mbox{\boldmath $x$}}_n^{(\ell)}$. To achieve this we propose the following choice: $$\begin{aligned}
  W_{n+1}=L_{n+1}^{(\ell)}\ {\textrm}{such that}\ L_{n+1}^{(\ell)}T_{n+1,n}^{(\ell)}=
\pmatrix{B_n^{(\ell)} \cr \mathbf{0}^T},   \label{bidiagonal_weight}\end{aligned}$$ where $B_n^{(\ell)}$ is an $n$-by-$n$ upper bidiagonal matrix of the form $$\begin{aligned}
B_{n}^{(\ell)}:=
\pmatrix{
t_{1,1}^{(\ell)}& t_{1,2}^{(\ell)} &        &       \cr
                & t_{2,2}^{(\ell)} & \ddots &       \cr
                &          &         \ddots &   t_{n-1,n}^{(\ell)} \cr
                &          &                &   t_{n,n}^{(\ell)}},\end{aligned}$$ and $L^{(\ell)}_{n+1}$ is lower triangular and will be specified below.

Next we derive recurrence formulas for updating the approximate solutions ${\mbox{\boldmath $x$}}_n^{(\ell)}$. From (\[WLSP\]) with the choice $W_{n+1}=L_{n+1}^{(\ell)}$ of (\[bidiagonal\_weight\]) it follows that $$\begin{aligned}
{\mbox{\boldmath $y$}}_n^{(\ell)}  &=& \arg\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\| g_1 {\mbox{\boldmath $e$}}_1- T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\Bigr\|_{(L_{n+1}^{(\ell)})^{\textrm{\scriptsize{H}}}L_{n+1}^{(\ell)}} \label{WLSP_B}\\
                &=& \arg\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \Bigl\| g_1  L_{n+1}^{(\ell)}{\mbox{\boldmath $e$}}_1-  L_{n+1}^{(\ell)}T_{n+1,n}^{(\ell)}{\mbox{\boldmath $z$}}_n^{(\ell)}\Bigr\|_{2}      \nonumber\\
                &=& \arg\min_{{\mbox{\boldmath $z$}}_n^{(\ell)}\in \mathbf{C}^n}\ \biggl\| 
\pmatrix{\widetilde{{\mbox{\boldmath $g$}}}_n^{(\ell)} \cr \widetilde{g}_{n+1}^{(\ell)}}
- \pmatrix{B_n^{(\ell)} \cr \mathbf{0}^{\textrm{\scriptsize{T}}}}
  {\mbox{\boldmath $z$}}_n^{(\ell)}\biggr\|_2, \ \pmatrix{\widetilde{{\mbox{\boldmath $g$}}}_n^{(\ell)} \cr \widetilde{g}_{n+1}^{(\ell)}}:=g_1 L_{n+1}^{(\ell)}{\mbox{\boldmath $e$}}_1. \nonumber\end{aligned}$$ From the above least squares problems we readily see that ${\mbox{\boldmath $y$}}_n^{(\ell)}=\ \bigr(B_n^{(\ell)}\bigl)^{-1}\widetilde{{\mbox{\boldmath $g$}}}_n^{(\ell)}$. Hence it follows from (\[sol\]) and using $(\widetilde{{\mbox{\boldmath $p$}}}_1\ \widetilde{{\mbox{\boldmath $p$}}}_2\ \cdots\ \widetilde{{\mbox{\boldmath $p$}}}_n):=V_n(B_n^{(\ell)})^{-1}$ that we have the following coupled two-term recurrence relations: $$\begin{aligned}
 \widetilde{{\mbox{\boldmath $p$}}}_{n}^{(\ell)}&=& ({\mbox{\boldmath $v$}}_{n}-t_{n-1,n}^{(\ell)}\widetilde{{\mbox{\boldmath $p$}}}_{n-1}^{(\ell)})/t_{n,n}^{(\ell)},   \label{x_rec1_B}\\
 {\mbox{\boldmath $x$}}_{n}^{(\ell)}&=&{\mbox{\boldmath $x$}}_{n-1}^{(\ell)}+\widetilde{g}_n^{(\ell)}\widetilde{{\mbox{\boldmath $p$}}}_n^{(\ell)}. \label{x_rec2_B}\end{aligned}$$ The cost per iteration is now $5Nm$. Substituting ${\mbox{\boldmath $p$}}_i^{(\ell)}=t_{i,i}^{(\ell)}\widetilde{{\mbox{\boldmath $p$}}}_i^{(\ell)}$ into the above recurrences, we have the even more efficient recurrence formulas. $$\begin{aligned}
 {\mbox{\boldmath $p$}}_{n}^{(\ell)}&=& {\mbox{\boldmath $v$}}_{n}-(t_{n-1,n}^{(\ell)}/t_{n-1,n-1}^{(\ell)}){\mbox{\boldmath $p$}}_{n-1}^{(\ell)}, \\
 {\mbox{\boldmath $x$}}_{n}^{(\ell)}&=&{\mbox{\boldmath $x$}}_{n-1}^{(\ell)}+(\widetilde{g}_n^{(\ell)}/t_{n,n}^{(\ell)}){\mbox{\boldmath $p$}}_n^{(\ell)}.\end{aligned}$$ By this reformulation, the cost becomes $4Nm+2m$.

Next, let us consider a choice for $L_{n+1}^{(\ell)}$ satisfying (\[bidiagonal\_weight\]). Let $F_{n}^{(\ell)}$ be an $n$-by-$n$ matrix for the form $$\begin{aligned}
  F_n^{(\ell)}(i):= 
  \pmatrix{
    I_{i-1}  &           &     & \cr
             &      1    &     & \cr
             &  f_{i}^{(\ell)}&   1 & \cr
             &           &     & I_{n-i-1}        
  },    \label{F}\end{aligned}$$ and let $T$ be tridiagonal. Then, by determining $f_i^{(\ell)}$ such that the $(i+1,i)$ entry of the matrix $F_n^{(\ell)}(i)T$ is zero, we can fulfill (\[bidiagonal\_weight\]) in the following way: $$\begin{aligned}
  F_{n+1}^{(\ell)}(n)F_{n+1}^{(\ell)}(n-1)\cdots F_{n+1}^{(\ell)}(1)T_{n+1,n}^{(\ell)}=
\pmatrix{
         B_n^{(\ell)} \cr
         \mathbf{0}^{\textrm{\scriptsize{T}}}\cr 
         }  \label{FN}.      \end{aligned}$$ From the above we see that $F_{n+1}^{(\ell)}(n)F_{n+1}^{(\ell)}(n-1)\cdots F_{n+1}^{(\ell)}(1)=L_{n+1}^{(\ell)}$, and thus $L_{n+1}^{(\ell)}$ and $L_{n}^{(\ell)}$ are related by $$\begin{aligned}
  L_{n+1}^{(\ell)} = F_{n+1}^{(\ell)}(n)\pmatrix{
         L_n^{(\ell)} & \mathbf{0}\cr  
         \mathbf{0}^{\textrm{\scriptsize{T}}}& 1\cr 
         } \    {\textrm}{for}\  n=2,3,\ldots,  \label{Relation_L}\end{aligned}$$ where $L_2^{(\ell)}=F_2^{(\ell)}(1).$ From the above we readily see that the matrices $L_{n+1}^{(\ell)}$ are lower triangular with all 1’s on the diagonals, and this property will be used for the proof of proposition 3.1 given later. The shifted QMR\_SYM method with the weight $(L_{n+1}^{(\ell)})^{\textrm{\scriptsize{H}}}L_{n+1}^{(\ell)}$ is referred to as shifted QMR\_SYM($B$). Its algorithm is given next.\
\
\
$$\begin{aligned}
&&{\hspace{-7mm}}\\[-15mm]
&&{\hspace{-7mm}}{\mbox{\boldmath $x$}}_{0}^{(\ell)}={\mbox{\boldmath $p$}}_{-1}^{(\ell)}=
      {\mbox{\boldmath $p$}}_{0}^{(\ell)}=\mathbf{0}, \ 
      {\mbox{\boldmath $v$}}_{1}={\mbox{\boldmath $b$}}/({\mbox{\boldmath $b$}}^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}})^{1/2},\ \widetilde{g}_1^{(\ell)}=({\mbox{\boldmath $b$}}^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}})^{1/2},\\
&&{\hspace{-7mm}}{\textrm}{{\bf for}}\ n=1,2, \ldots\ {\textrm}{{\bf do:}}\\
&&{\hspace{-7mm}}\quad {\textrm}{(The complex symmetric Lanczos process)}\nonumber\\
&&{\hspace{-7mm}}\quad \alpha_n ={\mbox{\boldmath $v$}}_n^{\textrm{\scriptsize{T}}}A{\mbox{\boldmath $v$}}_n,                           \nonumber \\
&&{\hspace{-7mm}}\quad \widetilde{{\mbox{\boldmath $v$}}}_{n+1}=
        A{\mbox{\boldmath $v$}}_n-\alpha_n{\mbox{\boldmath $v$}}_n-\beta_{n-1}{\mbox{\boldmath $v$}}_{n-1},              \nonumber \\
&&{\hspace{-7mm}}\quad \beta_{n} =(\widetilde{{\mbox{\boldmath $v$}}}_{n+1}^{\textrm{\scriptsize{T}}}\widetilde{{\mbox{\boldmath $v$}}}_{n+1})^{1/2},                  \nonumber \\
&&{\hspace{-7mm}}\quad {\mbox{\boldmath $v$}}_{n+1} = \widetilde{{\mbox{\boldmath $v$}}}_{n+1}/\beta_{n},          \nonumber \\
&&{\hspace{-7mm}}\quad t_{n-1,n}^{(\ell)}=\beta_{n-1},\ t_{n,n}^{(\ell)} =\alpha_{n}+\sigma_{\ell},\  t_{n+1,n}^{(\ell)}=\beta_{n},         \nonumber \\
&&{\hspace{-7mm}}\quad {\textrm}{(Solve weighted least squares problems)}                            \nonumber \\
&&{\hspace{-7mm}}\quad  {\textrm}{{\bf for}}\ \ell=1,2, \ldots, m\ {\textrm}{{\bf do:}}\\
&&{\hspace{-7mm}}\qquad {\textrm}{{\bf if}}\ \|{\mbox{\boldmath $r$}}_{n}^{(\ell)}\|_2/\|{\mbox{\boldmath $b$}}\|_2> \epsilon,\
            {\textrm}{{\bf then}}\\
&&{\hspace{-7mm}}\qquad \quad {\textrm}{{\bf for}}\ i=\max\{1,n-1\}, \ldots, n-1\ 
                   {\textrm}{{\bf do:}}\\
&&{\hspace{-7mm}}\qquad \qquad t_{i+1,n}^{(\ell)} =  f_{i}^{(\ell)} t_{i,n}^{(\ell)} +   t_{i+1,n}^{(\ell)},   \nonumber\\
&&{\hspace{-7mm}}\qquad \quad {\textrm}{{\bf end}}             \nonumber\\
&&{\hspace{-7mm}}\qquad  \quad f_n^{(\ell)}=-\frac{t_{n+1,n}^{(\ell)}}{t_{n,n}^{(\ell)}},   \nonumber \\   
&&{\hspace{-7mm}}\qquad  \quad t_{n+1,n}^{(\ell)}=0,                        \nonumber \\
&&{\hspace{-7mm}}\qquad \quad  \widetilde{g}_{n+1}^{(\ell)} =f_{n}^{(\ell)}\widetilde{g}_{n}^{(\ell)},  \nonumber\\
&&{\hspace{-7mm}}\qquad \quad  {\textrm}{(Update approximate solutions ${\mbox{\boldmath $x$}}_n^{(\ell)}$)}\nonumber\\
&&{\hspace{-7mm}}\qquad \quad {\mbox{\boldmath $p$}}_{n}^{(\ell)}= {\mbox{\boldmath $v$}}_{n}-(t_{n-1,n}^{(\ell)}/t_{n-1,n-1}^{(\ell)}){\mbox{\boldmath $p$}}_{n-1}^{(\ell)}, \nonumber\\
&&{\hspace{-7mm}}\qquad \quad {\mbox{\boldmath $x$}}_{n}^{(\ell)}={\mbox{\boldmath $x$}}_{n-1}^{(\ell)}+(\widetilde{g}_n^{(\ell)}/t_{n,n}^{(\ell)}){\mbox{\boldmath $p$}}_n^{(\ell)},  \nonumber \\
&&{\hspace{-7mm}}\quad \quad {\textrm}{{\bf end if}}       \nonumber    \\
&&{\hspace{-7mm}}\quad {\textrm}{{\bf end}}                \nonumber    \\&&{\hspace{-7mm}}\qquad {\textrm}{{\bf if}}\ \|{\mbox{\boldmath $r$}}_{n}^{(\ell)}\|_2 / \|{\mbox{\boldmath $b$}}\|_2\le \epsilon \ 
             {\textrm}{for all}\ \ell,\ {\textrm}{{\bf then} exit}.\\
&&{\hspace{-7mm}}{\textrm}{{\bf end}}                      \nonumber             \end{aligned}$$ Similar to Proposition 2.1, there is an efficient way to evaluate residual 2-norms as follows:

The $n$th residual 2-norms of the approximate solutions ${\mbox{\boldmath $x$}}_n^{(\ell)}$ for the shifted QMR\_SYM($B$) method are given by $$\begin{aligned}
 \|{\mbox{\boldmath $r$}}_n^{(\ell)}\|_2=|\widetilde{g}_{n+1}^{(\ell)}|\cdot\|{\mbox{\boldmath $v$}}_{n+1}\|_2\quad {\textrm}{for} \quad   \ell=1,2,\ldots,m.\end{aligned}$$

[*Proof*]{}. The proof is similar to that of Proposition 2.1. It follows from (\[QMR\_residual\]), (\[bidiagonal\_weight\]), (\[WLSP\_B\]), and recalling ${\mbox{\boldmath $y$}}_n^{(\ell)}=\ \bigr(B_n^{(\ell)}\bigl)^{-1}\widetilde{{\mbox{\boldmath $g$}}}_n^{(\ell)}$ that we have $$\begin{aligned}
{\mbox{\boldmath $r$}}_{n}^{(\ell)}=\widetilde{g}_{n+1}^{(\ell)}V_{n+1}\bigr(L_{n+1}^{(\ell)}\bigl)^{-1}{\mbox{\boldmath $e$}}_{n+1}. \end{aligned}$$ From (\[F\]) and (\[FN\]) $L_{n+1}^{(\ell)}$ is a lower triangular matrix with all 1’s on the diagonals. Thus, $\bigr(L_{n+1}^{(\ell)}\bigl)^{-1}$ is also a lower triangular matrix with all 1’s on the diagonals. It follows that $\bigr(L_{n+1}^{(\ell)}\bigl)^{-1}{\mbox{\boldmath $e$}}_{n+1}={\mbox{\boldmath $e$}}_{n+1}$, and thus we have ${\mbox{\boldmath $r$}}_{n}^{(\ell)}=\widetilde{g}_{n+1}^{(\ell)}V_{n+1}{\mbox{\boldmath $e$}}_{n+1}=\widetilde{g}_{n+1}^{(\ell)}{\mbox{\boldmath $v$}}_{n+1}$, which concludes the proof.

When we solve a large number of shifted systems, the computational cost of the residual 2-norms for the shifted QMR\_SYM($B$) method is much cheaper than the that for the shifted QMR\_SYM method since the former cost is of order $N$ and the latter is of order $Nm$ per iteration step.

Observing algorithms of the shifted QMR\_SYM and the shifted QMR\_SYM($B$) methods, we see that the work for the weighted least squares problems and for updating approximate solutions in the shifted QMR\_SYM($B$) method is lower than that in the shifted QMR\_SYM method. For the case $m\gg 1$, the most time-consuming part is updating the approximate solutions. In this case, the shifted QMR\_SYM($B$) method may be more efficient than the shifted QMR\_SYM method since the shifted QMR\_SYM($B$) method requires $4Nm+2m$ operations per iteration step for the update while the shifted QMR\_SYM method needs $6Nm+3m$ operations. This difference will be clearer when we use the results of Propositions 2.1 and 3.1 as stopping criterion. On the other hand, in terms of number of iterations, under a certain assumption, the convergence of the shifted QMR\_SYM($B$) method is worse than that of the shifted QMR\_SYM method but not worse than that of the shifted COCG method, as described by the following proposition:

Under the assumption of Proposition 2.2, the shifted QMR\_SY\
M($B$) method (Algorithm 3.2) enjoys the following properties:

-   all matrix-vector multiplications can be done in real arithmetic;

-   if breakdown does not occur and each matrix $T_n+\sigma_\ell I_n$is nonsingular, then $\bigr\|{\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SQ($B$)}}}\bigl\|_2=\bigr\|{\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}\bigl\|_2\ge \bigr\|{\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SQ}}}\bigl\|_2 \ {\textrm}{for}\   \ell=1,2,\ldots m$, where the superscripts SQ($B$), SCOCG, and SQ are short for shifted QMR\_SYM($B$), shifted COCG, and shifted QMR\_SYM, respectively;

-   $ \|{\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SQ($B$)}}}\|_2=|\widetilde{g}_{n+1}^{(\ell)}|\  {\textrm}{for} \  \ell=1,2,\ldots,m, \  n\ge 0.$

[*Proof*]{}. The proof of ($I$) is the same as that of Proposition 2.2, and is based on the fact that under the assumption the complex symmetric Lanczos process generates real basis vectors.

Next, we give a proof of ($I\hspace{-0.5mm}I$). The $n$th residuals of the shifted COCG method [@Shifted_COCG] belong to ${\mbox{\boldmath $b$}}-(A+\sigma_\ell I)K_n(A+\sigma_\ell I,{\mbox{\boldmath $b$}})$. Hence, each ${\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}$ can be written as ${\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}={\mbox{\boldmath $b$}}-(A+\sigma_\ell I)V_n{\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}$, where $V_n$ is the same matrix as in (\[comp\_sym\_Lanczos\_Mat\]) Since each ${\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}$ is orthogonal to each subspace $K_n(\overline{A}+\overline{\sigma}_\ell I,\overline{{\mbox{\boldmath $b$}}})${=$K_n(\overline{A},\overline{{\mbox{\boldmath $b$}}})\}$, i.e., ${\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}} \perp K_n(\overline{A},\overline{{\mbox{\boldmath $b$}}})$, we have $V_n^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}}-V_n^{\textrm{\scriptsize{T}}}(A+\sigma_\ell I)V_n{\mbox{\boldmath $y$}}^{(\ell), {\textrm}{\scriptsize{SCOCG}}}=0$, and thus it follows from (\[comp\_sym\_Lanczos\_Mat\]) that we have the relation ${\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}} = \{V_n^{\textrm{\scriptsize{T}}}(A+\sigma_\ell I)V_n\}^{-1}V_n^{\textrm{\scriptsize{T}}}{\mbox{\boldmath $b$}}= g_1(T_n^{(\ell)})^{-1}{\mbox{\boldmath $e$}}_1$, where $T_n^{(\ell)}:=T_n+\sigma_\ell I_n$. Since the shifted QMR\_SYM($B$) method has the form (\[QMR\_residual\]), it is sufficient to show that ${\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}={\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SQ($B$)}}}$. From (\[WLSP\_B\]) and (\[Relation\_L\]) it follows that ${\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SQ($B$)}}}=(B_n^{(\ell)})^{-1}\widetilde{{\mbox{\boldmath $g$}}}_n^{(\ell)}=g_1(B_n^{(\ell)})^{-1}[I_n |\mathbf{0}]L_{n+1}^{(\ell)}[{\mbox{\boldmath $e$}}_1^{\textrm{\scriptsize{T}}}| 0]^{\textrm{\scriptsize{T}}}=g_1(B_n^{(\ell)})^{-1}L_{n}^{(\ell)}{\mbox{\boldmath $e$}}_1=g_1\{(L_{n}^{(\ell)})^{-1}B_n^{(\ell)}\}^{-1} {\mbox{\boldmath $e$}}_1.$ Since from (\[FN\]) and (\[Relation\_L\]) we can readily confirm the relation $L_{n}^{(\ell)}T_n^{(\ell)}=B_n^{(\ell)}$, we have ${\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SQ($B$)}}}=g_1(T_n^{(\ell)})^{-1} {\mbox{\boldmath $e$}}_1,$ which is the same as ${\mbox{\boldmath $y$}}_n^{(\ell), {\textrm}{\scriptsize{SCOCG}}}.$ The inequality in ($I\hspace{-0.5mm}I$) follows from Proposition 2.2 since under the given assumption the shifted QMR\_SYM method enjoys the minimal residual property.

Finally, we give a proof of ($I\hspace{-0.7mm}I\hspace{-0.7mm}I$). If follows from the proof of ($I$) that $\|{\mbox{\boldmath $v$}}_i\|_2=1$ for all $i$. Thus from Proposition 3.1 we have $ \|{\mbox{\boldmath $r$}}_n^{(\ell), {\textrm}{\scriptsize{SQ$(B)$}}}\|_2=|\widetilde{g}_{n+1}^{(\ell)}|\cdot\|{\mbox{\boldmath $v$}}_{n+1}\|_2=|\widetilde{g}_{n+1}^{(\ell)}|\ {\textrm}{for} \   \ell=1,2,\ldots,m,\ n\ge 0.$

In the property ($I\hspace{-0.7mm}I$) of Proposition 3.2, breakdown may occur due to the choice (\[F\]) of the weighted least squares problems.

From proposition 3.2 we see that in terms of the number of iteration steps the shifted QMR\_SYM($B$) method never converges faster than the shifted QMR\_SYM method but it converges at the same iteration step as the shifted COCG method does. Since the efficiency of the shifted COCG method has been shown already and for the case $m\gg 1$ the computational cost of the shifted QMR\_SYM$(B)$ method is much less than that of the shifted QMR\_SYM method, the shifted QMR\_SYM$(B)$ method can also be useful. This will be supported by some numerical examples in the next section.

Numerical examples
==================

In this section, we report on some numerical examples with the shifted COCG method, the shifted QMR\_SYM method (Algorithm 2.1), and the shifted QMR\_SYM($B$) method (Algorithm 3.2). We evaluate these methods in terms of computational time. All tests were performed on a workstation with a 2.6GHz AMD Opteron(tm) processor 252 using double precision arithmetic. Codes were written in Fortran 77 and compiled with g77 -$O$3. In all cases the stopping criterion was set as $\epsilon=10^{-12}$.\
\
3.1. Example 1\
\
The first problem comes from the electronic structure calculation of a bulk Si(001) with 512 atoms in [@Shifted_COCG], which is written as follows: $$\begin{aligned}
(\sigma_\ell I - H){\mbox{\boldmath $x$}}^{(\ell)} = {\mbox{\boldmath $e$}}_1, \ \  \ell=1,2, \ldots,m,\end{aligned}$$ where $\sigma_\ell=0.400 + (\ell-1+i)/1000$, $H\in R^{2048\times 2048}$ is a symmetric matrix with 139264 entries, ${\mbox{\boldmath $e$}}_1=(1,0,\ldots,0)^{\textrm{\scriptsize{T}}}$, and $m=1001$. Since the shifted COCG method requires the choice of a seed system, we have chosen the optimal seed ($\ell=714$) in this problem; otherwise some linear systems will remain unsolved.

Figure 4.1 shows the number of iterations of each method to solve the $\ell$th shifted linear systems. For example, in Fig. 4.1, the number of iterations for the shifted COCG method at $\ell=600$ is 150, which means the shifted COCG method required 150 iterations to obtain the (approximate) solutions of the 600th shifted linear system, i.e., $(\sigma_{600} I - H){\mbox{\boldmath $x$}}^{(600)} = {\mbox{\boldmath $e$}}_1$.

![Number of iterations for the shifted COCG method, the shifted QMR\_SYM method, and the shifted QMR\_SYM($B$) method versus the serial number of the shifted linear systems. ](iteration.eps){width="0.8\linewidth" height="0.5\linewidth"}

![[Log]{}$_{10}$ of the relative residual 2-norms versus the number of iterations of the shifted COCG method, the shifted QMR\_SYM method, and the shifted QMR\_SYM($B$) method for the shifted linear system with $\ell=701, i.e., \sigma_{701}=1.100+0.001i$.](history.eps){width="0.8\linewidth" height="0.5\linewidth"}

From Fig. 4.1 we make three observations: first, the three methods required almost the same number of iterations at each $\ell$; second, in terms of number of iterations, the shifted QMR\_SYM method often converged slightly faster than the other two methods. This phenomenon is closely related to Proposition 2.2, which will become clearer later; third, for each method the required number of iterations depends highly on the shift parameters $\sigma_\ell$. This result may come from varying eigenvalues of the coefficient matrices $\sigma_\ell I -H$ since if we choose $\sigma_\ell$ close to an eigenvalue of $H$, then $\sigma_\ell I -H$ is close to singular. Conversely, from the shape of Fig. 4.1 one may obtain a partial distribution of eigenvalues of $H$.

One of the residual 2-norm histories for the three methods is given in Fig. 4.2. From Fig. 4.2 we see that the relative residual 2-norm of the shifted QMR\_SYM method decreases monotonically and at every iteration step the norm is less than those of the other two methods. Hence we can say that the property $(I\hspace{-0.5mm}I)$ of Proposition 2.2 is experimentally supported by this history. We make another observation in the histories of the shifted COCG method and the shifted QMR\_SYM($B$) method: During the first about fifty iterations, the two methods show the same histories. After that their histories varies gradually. Hence we see that the property $(I\hspace{-0.5mm}I)$ of Proposition 3.2 is also experimentally supported by these histories.

![Required CPU time given in seconds versus the number of shifted linear systems for each iterative method.](Si512_CPU_time.eps){width="0.8\linewidth" height="0.5\linewidth"}

![The ratio of each computational time to the one of the shifted COCG method versus the number of shifted linear systems.](Si512_CPU_time_ratio.eps){width="0.8\linewidth" height="0.5\linewidth"}

Each computational time of the three methods is given in Fig. 4.3, where the horizontal axis denotes the number $m$ of shifted linear systems that are solved. For example, in Fig. 4.3, the computational time of the shifted COCG method at $m=200$ is about 0.76 \[sec.\], which means that it required about 0.76 \[sec.\] to solve the shifted linear systems: $\{(0.400+0.001i) I - H\}{\mbox{\boldmath $x$}}^{(1)} = {\mbox{\boldmath $e$}}_1, \{(0.401+0.001i) I - H\}{\mbox{\boldmath $x$}}^{(2)} = {\mbox{\boldmath $e$}}_1, \ldots, \{(0.599+0.001i) I - H\}{\mbox{\boldmath $x$}}^{(200)} = {\mbox{\boldmath $e$}}_1$. From Fig. 4.3 we see that as the number $m$ grows larger, the shifted QMR\_SYM method required more CPU time than the other two methods. On the other hand, the shifted QMR\_SYM($B$) method required almost the same CPU time as the shifted COCG method. This phenomena can be attributed to the computational costs of updating approximate solutions for each method since there are three facts: first, we know from Fig. 4.1 that the three methods required almost the same iterations; second, the shifted QMR\_SYM($B$) method requires almost the same computational costs as the shifted COCG method, while the shifted QMR\_SYM method tends to require more computational costs per iteration than the other two methods; third, for large $m$, updating approximate solutions is one of the most time-consuming parts. The second and the third facts have been discussed already in the previous two sections.

In Fig. 4.3 we can see little about the properties of the three methods for small $\ell$. We therefore show the ratio of each computational time to the computational time of the shifted COCG method in Fig. 4.4. We see from Fig. 4.4 that in terms of the ratio of CPU times the shifted QMR\_SYM method and the shifted QMR\_SYM($B$) method converged much faster than the shifted COCG method when the number of shifted linear systems is small, say, $m<200$. This can be explained in the following way: for small $m$, updating approximate solutions does not affect the CPU time so much. Other operations such as matrix-vector multiplications are now the most time-consuming parts since the three methods required almost the same number of iterations, see Fig. 4.1. From Proposition 2.2 $(I)$ and Proposition 3.2 $(I)$ we know that in this case the shifted QMR\_SYM method and the shifted QMR\_SYM($B$) method require real matrix-real vector multiplications. On the other hand, the shifted COCG method requires real matrix-complex vector multiplications. Moreover, dot products and vector additions of the complex symmetric Lanczos process used in the shifted QMR\_SYM method and the shifted QMR\_SYM($B$) method can be done in real arithmetic. Hence, the two methods converged much faster than the shifted COCG method.\
\
3.2. Example 2\
\
The second problem comes from the electronic structure calculation of bulk fcc Cu with 1568 atoms in [@Shifted_COCG] and is given as follows: $$\begin{aligned}
(\sigma_\ell I - H){\mbox{\boldmath $x$}}^{(\ell)} = {\mbox{\boldmath $e$}}_1, \ \  \ell=1,2,\ldots,m,\end{aligned}$$ where $\sigma_\ell=-0.5+(\ell-1+i)/1000$, $H\in R^{14112\times 14112}$ is a symmetric matrix with 3924704 entries, ${\mbox{\boldmath $e$}}_1=(1,0,\ldots,0)^{\textrm{\scriptsize{T}}}$, and $m=1501$.

![The ratio of each computational time to the one of the shifted COCG method versus the number of shifted linear systems.](Cu1568_CPU_time.eps){width="0.8\linewidth" height="0.5\linewidth"}

![The ratio of each computational time to the one of the shifted COCG method versus the number of shifted linear systems.](Cu1568_CPU_time_ratio.eps){width="0.8\linewidth" height="0.5\linewidth"}

Each computational time of the three methods for solving the $m$ shifted linear systems is shown in Fig. 4.5, where horizontal and vertical axes are the same as those used in Fig 4.3. The ratio of each computational time to the computational time of the shifted COCG method is shown in Fig. 4.6. From Figs. 4.5, 4.6 we see that the size of this matrix is about 7 times larger than before, and the three methods show similar tendencies to those seen in the previous example.

Concluding remarks
==================

In this paper, the shifted QMR\_SYM method was described as a specialization of the QMR method for general non-Hermitian shifted linear systems [@Freund]. The method has an advantage over the shifted COCG method in that there is no need to choose a suitable seed system. On the other hand, we have found that for a large number of shifted linear systems, the most time-consuming part of the shifted QMR\_SYM method is updating approximate solutions, and this cost is higher than that of the shifted COCG method. We therefore have proposed the weighted quasi-minimal residual approach with a specific weight for reducing the computational cost for updating approximate solutions. The resulting method, shifted QMR\_SYM($B$), also does not require to choose a suitable seed system, which is an advantage over the shifted COCG method. From numerical experiments we have learned that shifted QMR\_SYM and QMR\_SYM($B$) are competitive in comparison to the shifted COCG method. In particular, QMR\_SYM($B$) can be the method of choice for solving complex symmetric shifted linear systems with a large number of shifts that arise from large-scale electronic structure theory. In future work, numerical tests for general complex symmetric shifted linear systems will be done to evaluate the performance of the method.\

Acknowledgments {#acknowledgments .unnumbered}
===============

We wish to express our gratitude to Roland W. Freund (UC Davis) for his fruitful comments during the conference at Harrachov in 2007. We are grateful to an anonymous referee for useful comments that substantially enhanced the quality of this paper. Finally, we would like to thank Martin H. Gutknecht (ETH Zurich) for careful reading of the manuscript.\
\

[10]{}

, [*Necessary and sufficient conditions for the existence of a conjugate gradient method*]{}, SIAM J. Numer. Anal., 21 (1984), pp. 352-362.

, [*On conjugate gradient type methods and polynomial preconditioners for a class of complex non-Hermitian matrices*]{}, Numer. Math., 57 (1990), pp. 285-312.

, [*Conjugate gradient-type methods for linear systems with complex symmetric coefficient matrices*]{}, SIAM J. Sci. Stat. Comput., 13 (1992), pp. 425-448.

, [*Solution of Shifted Linear Systems by Quasi-Minimal Residual Iterations*]{}, Numerical Linear Algebra, L. Reichel, A. Ruttan and R. S. Varga eds., W. de Gruyter, (1993), pp. 101-121.

, [*QMR: A quasi-minimal residual method for non-Hermitian linear systems*]{}, Numer. Math., 60 (1991), pp. 315-339.

, [*BiCGStab($\ell$) for families of shifted linear systems*]{}, Computing 70 (2003), pp. 87-109.

, [*Restarted GMRES for shifted linear systems*]{}, SIAM J. Sci. Comput. 19 (1998), pp. 15-26.

, [*Matrix computations*]{}, 3rd ed., The Johns Hopkins University Press, Baltimore and London, 1996.

, [*An iteration method for the solution of the eigenvalue problem of linear differential and integral operators*]{}, J. Res. Nut. Bur. Standards, 45 (1950), pp. 255-282.

, [*Restarted full orthogonalization method for shifted linear systems*]{}, BIT Numerical Mathematics, 43 (2003), pp. 459-466.

, [*Recent computational developments in Krylov subspace methods for linear systems*]{}, Num. Lin. Alg. with Appl., 14 (2007), pp. 1-59.

, [*A projection method for generalized eigenvalue problems using numerical integration*]{}, J. Comput. Appl. Math., 159 (2003), pp. 119-128.

, [*A COCR method for solving complex symmetric linear systems*]{}, J. Comput. Appl. Math., 199 (2007), pp. 297-303.

, [*A numerical method for calculating the Green’s function arising from electronic structure theory*]{}, in: Frontiers of Computational Science, eds. Y. Kaneda, H. Kawamura and M. Sasai, Springer-Verlag, Berlin/Heidelberg, 2007, pp. 189-195.

, [*Linear algebraic calculation of Green’s function for large-scale electronic structure theory*]{}, Phys. Rev. B 73, 165108 (2006), pp. 1-9.

, [*A Petrov-Galerkin type method for solving $A{\mbox{\boldmath $x$}}={\mbox{\boldmath $b$}}$, where $A$ is symmetric complex*]{}, IEEE Trans. Mag. 26(2) (1990), pp. 706-708.

[^1]: Department of Computational Science and Engineering, Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan ([{sogabe,zhang}@na.cse.nagoya-u.ac.jp]{}).

[^2]: \(1) Department of Applied Mathematics and Physics, Tottori University, Japan (2) Core Research for Evolutional Science and Technology, Japan Science and Technology Agency (CREST-JST), 4-1-8 Honcho, Kawaguchi-shi, Saitama 332-0012, Japan ([hoshi@damp.tottori-u.ac.jp]{}).

[^3]: \(1) Center for Research and Development of Higher Education, The University of Tokyo, Hongo, 7-3-1, Bunkyo-ku, Tokyo 113-8656, Japan (2) Core Research for Evolutional Science and Technology, Japan Science and Technology Agency (CREST-JST), 4-1-8 Honcho, Kawaguchi-shi, Saitama 332-0012, Japan ([fujiwara@coral.t.u-tokyo.ac.jp]{}).

[^4]: This work was partially supported by KAKENHI (Grant No. 18760063, 19560065).
