# A Survey of Low-Rank Updates of Preconditioners for Sequences of Symmetric Linear Systems
## 7 Numerical Results
### Fe Discretization of a Parabolic Equation  
Consider the parabolic PDE  
\[S_{h}\frac{\partial u(\vec{x},t)}{\partial t} = \
abla\left(K(\vec{x})\
abla u(\vec{x},t)\right)+\text{ BCs}.\]  
\(u\) is the pressure of the fluid, \(S_{h}\) the elastic storage, \(K\) the (SPD) hydraulic conductivity tensor. This equation is used to model water pumping in the Beijing plan to predict land subsidence in that area [45]. The 3D computational domain turns out to be very heterogeneous and geometrically irregular as sketched in Figure 3:  
FE discretization in space needs highly irregular tetrahedra, which, combined with the high heterogeneity of the medium, gives raise to an ill-conditioned discretized diffusion matrix. After FE discretization a system of ODEs is left:  
\[B\frac{\partial\mathbf{u}(t)}{\partial t} = -A\mathbf{u}(t)+\mathbf{b}(t).\]  
where \(A\) is the SPD stiffness matrix, \(B\) is the SPD mass matrix. The right-hand-side account for source terms and Neumann boundary conditions. Using a time marching scheme (e.g., Crank-Nicolson for stiff problems) yields a sequence of linear systems  
\[(A+\sigma_{k}B)\mathbf{u}_{k}=\mathbf{b}_{k},\quad\sigma_{k}=\frac{2}{\Delta t _{k}},\qquad k=1,Nstep. \tag{12}\]  
Figure 3: 3D domain and triangulation. The discretized FE matrices have both \(n\) = 589,661 and a number of nonzeros \(nz\) = 8,833,795. We now report the results of a complete simulation which took \(Nstep=32\) steps. The near steady state is reached after \(T\) = 36,000 days. Initially we set \(\Delta t_{1}=60\) (days). Then we incremented the timesteps as  
\[\Delta t_{k}=\min(1.2\Delta t_{k-1},7300,T-\Delta t_{k-1})\Longrightarrow\ \sigma_{k}=\frac{2}{\Delta t_{k}},\qquad k=2,\ldots,Nstep.\]  
This formula allows the timestep size to constantly increase and to take its largest value approaching the steady state solution. The initial preconditioner \(P_{0}\) was evaluated employing the Matlab function ichol with  
1. \(\texttt{droptol}=10^{-5}\) applied to \(A_{1}=A+5B\) to enhance diagonal dominance; this resulted in a triangular Cholesky factor with density 5.18 and
2. \(\texttt{droptol}=10^{-4}\) applied to \(A_{2}=A+20B\) (density 2.83). We solved the first linear system to a high accuracy to assess as accurately as possible the \(p=10\) leftmost eigenvectors. Using the residual \(\delta_{j}=\frac{\|P_{0}A\mathbf{w}_{j}-\mu_{j}\mathbf{w}_{j}\|}{\|\mathbf{w}_ {j}\|}\) to measure the accuracy of the approximate eigenvectors we found:  
* Case (1). 179 PCG iterations for the first linear system; \(\delta_{j}\in[1.1\,10^{-3},5.8\,10^{-3}]\). * Case (2). 370 PCG iterations for the first linear system; \(\delta_{j}\in[1.1\,10^{-3},3.3\,10^{-3}]\). This set of vectors is then used to construct deflation, tuned and spectral preconditioners for all the subsequent linear systems in the sequence (with exit test on relative residual and tolerance tol\(=1e-10\)). The results are summarized in Table 6. The eigenvectors estimated during the first linear system are adequate to accelerate the subsequent iterations for both the tuned and spectral approaches which provide an improvement in both iteration number and CPU time. Note also that the time per iteration is only slightly increased with respect to the \"no update\" case. Instead, the PCG residuals obtained with the deflation preconditioner in some instances stagnate before reaching the exit test. The lack of robustness of the deflation approach is also observed in Reference [5] where the authors suggest a reorthogonalization of the residuals at the price of increasing costs per iteration and also Reference [26] in the framework of iterative eigensolvers. Regarding tuned preconditioners, SR1 is definitely superior to the BFGS approach due to its lower computational cost.