&        &                \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, 