be used for any iterative solver generating approximations $\bxk$ and $\byk$, and residuals $\brk$ and $\bsk$, $k=1,2,\dots$. If the norms of the residual vectors $\brk$ and $\bsk$ tend to be small, one can expect that $$\begin{aligned}
  \label{BB5}
  | \bskT\mA^{-1} \brk | \ll |\bcT \mA^{-1} \brk| \qquad \mbox{and}\qquad
  |\bskT\mA^{-1} \brk| \ll |  \bskT\mA^{-1} \bb|.\end{aligned}$$ In other words, one can expect that the approximation is more accurate than since the corresponding error terms tend to be smaller.

Approximating the quantity of interest using BiCG iterates {#sec:bicg}
----------------------------------------------------------

As mentioned in Introduction, our aim is to solve systems by an iterative method, which allows to solve the primal and dual problems simultaneously. Since we intend to solve large and sparse systems, we need to pick a method with low memory requirements. Then, a natural choice is to use the [*preconditioned bi-conjugate gradient*]{} (BiCG) method (Algorithm \[alg:BiCG\]) introduced in [@Fl1976]; see also [@B:BaBeCh1994; @StTi2011]. Let $\mP$ be a suitable preconditioner (its choice depends on the particular discretization scheme). BiCG is a short-term recurrence Krylov subspace method which generates (if no breakdown occurs) approximations $\bxk \in \bx_0 + \mP^{-1}\mathcal{K}_k (\mA \mP^{-1},\br_0)$ and $\byk
\in \by_0 + \mP^{-\T}\mathcal{K}_k (\mA^{\T}\mP^{-\T},\bs_0)$ such that $$\label{eqn:bicg}
    \brk \perp \mP^{-\T}\mathcal{K}_k (\mA^{\T}\mP^{-\T},\bs_0),\qquad 