the error vector $$\boldsymbol{\delta}_{{\boldsymbol u},h_{k+1}}^{k}:=\nabla f(\boldsymbol{\rm u}^{k+1}_{h_{k+1}})+M_{h_{k+1}}^{T}\boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm \lambda}^{k}_{h_{k}}+\sigma M_{h_{k+1}}(\boldsymbol{\rm u}^{k+1}_{h_{k+1}}-\boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm z}^{k}_{h_{k}})$$      satisfies $\|\boldsymbol{\delta}_{{\boldsymbol u},h_{k+1}}^{k}\|_{2}\leq \xi_{k+1}$. Compute $\boldsymbol{\rm z}^{k}$ as follows: $$\boldsymbol{\rm z}^{k+1}_{h_{k+1}}={\rm argmin} \ g(\boldsymbol{\rm z})+\langle \boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm \lambda}^{k}_{h_{k}},M_{h_{k+1}}(\boldsymbol{\rm z}-\boldsymbol{\rm u}^{k+1}_{h_{k+1}})\rangle+\dfrac{\sigma}{2}\|\boldsymbol{\rm u}^{k+1}_{h_{k+1}}-\boldsymbol{\rm z}\|_{M_{h_{k+1}}}^{2}.$$ Compute $$\boldsymbol{\rm \lambda}^{k+1}_{h_{k+1}}=\boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm \lambda}^{k}_{h_{k}}+\tau\sigma(\boldsymbol{\rm u}^{k+1}_{h_{k+1}}-\boldsymbol{\rm z}^{k+1}_{h_{k+1}}).$$ If a termination criterion is met, stop; else, set $k:=k+1$ and go to Step 1.

Let $\boldsymbol{\rm y}^{k+1}:=K_{h}^{-1}M_{h}(\boldsymbol{\rm u}^{k+1}+\boldsymbol{\rm y}_{r}),\
\boldsymbol{\rm p}^{k+1}:=K_{h}^{-1}M_{h}(\boldsymbol{\rm y}_{d}-\boldsymbol{\rm y}^{k+1})$ denotes the discretized state and the discretized adjoint state respectively. Then the $\boldsymbol {\rm u}$-subproblem at the $k$-th iterate is equivalent to solving the following linear system: $$\left[                
  \begin{array}{ccc}  
    M_{h} & 0 & K_{h}\\  
    0 &  (\alpha+\sigma)M_{h} & -M_{h}\\ 
    K_{h} & -M_{h} & 0\\ 
  \end{array}
\right]
\begin{array}{ccc}  
\left[       
  \begin{array}{ccc}           
    \boldsymbol{\rm y}^{k+1} \\  
    \boldsymbol{\rm u}^{k+1}\\ 
    \boldsymbol{\rm p}^{k+1}\\ 
\end{array}
\right]
\end{array}
=
\left[ 
\begin{array}{ccc}           
    M_{h}\boldsymbol{\rm y}_{d} \\  
    M_{h}(-\boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm \lambda}^{k}+\sigma \boldsymbol{\rm I}_{h_{k+1}}{\boldsymbol{\rm z}}^{k})\\ 
    M_{h}\boldsymbol{\rm y}_{r}\\ 
\end{array}   
\right].$$ Since $\boldsymbol{\rm p}^{k+1}=(\alpha+\sigma)\boldsymbol{\rm u}^{k+1}+\boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm \lambda}^{k}-\sigma \boldsymbol{\rm I}_{h_{k+1}}\boldsymbol{\rm z}^{k}$, we eliminate the variable $\boldsymbol{\rm p}$, then the above problem can be rewritten in the following reduced form without any computational cost: $$\label{equ:linear}
\left[                
  \begin{array}{ccc}  
    M_{h} &  K_{h}(\alpha+\sigma) \\  
    -K_{h} &  M_{h} \\ 
  \end{array}
\right]
\begin{array}{ccc}  
\left[       
  \begin{array}{ccc}           
    \boldsymbol{\rm y}^{k+1} \\  
    \boldsymbol{\rm u}^{k+1}\\ 
\end{array}
\right]
\end{array}
=
\left[ 
\begin{array}{ccc}           
  M_{h}\boldsymbol{\rm 