- Generated by
1.8.13
|
HLIBpro
2.6
|
While the technical details of solving and preconditioning are discussed in Solvers, this section should comment on the practical issues related to the various preconditioners.
As a benchmark, we will use the same problem as in the Boundary Element Methods example, i.e., the Helmholtz SLP operator on the unit sphere. The wavenumber \(\kappa\) was chosen to be 16.
The source code for setting up the matrix and RHS is identical to the Boundary Element Methods example.
For this specific test, we use a grid with 62.466 vertices and 124.928 triangles. With a constant ansatz, this results in an H-matrix of dimension \(124.928 \times 124.928\).
Furthermore, the matrix was built unsymmetric, i.e., both lower and upper triangular part are constructed, to prevent algorithmic properties of symmetric matrices, e.g., during matrix-vector solves, to influence the numerical results.
For comparison with the timings below, we need the time to build the H-matrix, which is 75.9 seconds for a block-wise accuracy of \(10^{-4}\) with ACA+.
All tests below are performed using a single Intel Xeon E7-8857 CPU with 12 cores and 3GHz.
First, the various solvers in 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 should try to solve the system without any preconditioner. The stop criterion is chosen with a maximal number of 1000 steps and a reduction of the residual norm of \(10^{-6}\) compared to the initial residual norm. Furthermore, the lower bound for the absolute residual norm was set to \(10^{-6}\):
The results are
The fastest solver is TFQMR, which managed to reach the wanted residual norm within 575 steps. However, it needed 58.12 seconds to do so, which is almost as long as for building the matrix.
The next test will use the Jacobi preconditioner, first the point-wise version:
The results now look as
The GMRES solver now is able to reduce the residual within 193 steps in 22.97 seconds. Also improved are the convergence rates of BiCGStab and TFQMR. CGS now is not able to converge, resulting in absent data.
However, point-wise Jacobi only uses the diagonal entries of \(A\) as an approximation. In the next test, the diagonal blocks of \(A\) are used instead:
The third argument to the TJacobi constructor defines the block size and since no block smaller than diagonal blocks are supported by 𝖧𝖫𝖨𝖡𝗉𝗋𝗈, a 2 results in using all diagonal blocks for a block-wise Jacobi preconditioner.
Instead of a blocksize value of 2, one could also use \(n_{\min}\) if known, since this is equivalent as it defines the size of the smallest diagonal blocks.
The results are:
The convergence rate actually decreases for GMRES and BiCGStab. Only TFQMR benefits a little.
The reason for this behaviour is, that although going from point-wise diagonal to block-diagonal with blocksize \(n_{\min}\) will increase the number of coefficients in the preconditioner, but the change in the overall approximation of \(A^{-1}\) is still small, and hence the quality of the preconditioner does not change enbough to guarantee a better convergence for this example.
To improve this, we increase the diagonal blocksize to 100:
which gives the following results:
For comparison, also a blocksize of 500 is used:
Since the approximation to \(A\) and hence also to \(A^{-1}\) is getter better with an increased blocksize, the convergence rates are also decreasing.
However, an increased block size also has drawbacks:
The increased time for matrix-vector multiplications can be seen for the BiCGStab and the TFQMR solver, which have similar solve times for both block sizes, although the number of iteration steps is reduced.
For point-wise Jacobi and for a block-size equal to \(n_{\min}\), the setup was actually negligible. So lets correct the above data by the measured time for setting up the preconditioner. First for a blocksize of \(100\):
and also for a blocksize of 500:
The relatively large increase in the setup time is due to H-inversion for the diagonal blocks, as performed by the TJacobi class. For small blocks, this is comparable to dense inversion, which is very fast for small matrices. However, depending on the H-block structure, a more than linear increase may be visible for larger blocks (before H-complexity is dominant).
TJacobi also accepts a fourth argument defining the accuracy of the H-inversion, which defaults to exact inversion (up to machine precision):
In this example however, a relaxed accuracy only leads to very little change in the setup and solve times.
This is a typicall example that a more advanced preconditioner does not necessarily results in an overall performance increase. Only if the reduction in iteration steps is large enough to actually compensate the increase in the more costly setup phase, will the overall time also be reduced.
In addition to the diagonal of \(A\), the Gauss-Seidel preconditioner also uses the lower (forward GS), upper (backward GS) or both (symmetric GS) triangular part(s) of the matrix. This gives as the following preconditioners with \(A = D - E - F\) and \(E\) being the strict (block) lower triangular part of \(A\) and \(F\) the strict upper triangular part respectively:
Let us test all three versions of Gauss-Seidel for the Helmholtz model problem.
Forward-GS:
Backward-GS:
Symmetric-GS:
Only Symmetric-GS shows a significant improvement over block-wise Jacobi in terms of convergence rates, with less that 100 iteration steps to solve the Helmholtz problem.
However, the solve time is significantly higher compared to Jacobi which is the result of the matrix solve algorithm in 𝖧𝖫𝖨𝖡𝗉𝗋𝗈. For one, solving a lower/upper triangular H-matrix takes much longer than matrix-vector multiplication with a block-diagonal matrix since much more data is involved. Furthermore this matrix solve algorithm does not scale as efficiently with the number of CPUs as does matrix-vector multiplication, especially for a block-diagonal matrix.
Already Block-Jacobi and Gauss-Seidel made use of the block structure of H-matrices and the correponding H-arithmetic. Especially Block-Jacobi used H-matrix inversion for the diagonal blocks.
Using H-LU factorisation, a better approximation to \(A\) can be computed by using the whole matrix instead of just the (block-) diagonal, which then in theory should serve as an excellent preconditioner. In fact, an exact (up to machine precision) computation of \(A^{-1}\) would result in a direct solver, i.e., only a single matrix-vector multiplication (or matrix-vector solve) would be needed for solving the equation system.
H-arithmetic allows us to control the accuracy of the approximation, which thereby permits us to optimize the overall solution process, e.g., invest more time into setup (high accuracy H-LU) or more time in the iteration process (low accuracy H-LU) depending on our problem to solve.
For the Helmholtz problem, we will increase the block-wise H-arithmetic accuracy \(\varepsilon\) from \(1\) to \(10^{-4}\).
Except for the accuracy setting, the source for all examples is:
The numerical results are:
\(\varepsilon = 1\)
\(\varepsilon = 10^{-2}\)
\(\varepsilon = 10^{-4}\)
As expected, with a better H-accuracy, the convergence rate also becomes better, until only two iteration steps are needed for \(\varepsilon = 10^{-4}\). However, the setup time for computing the H-LU preconditioner also increases significantly with a smaller \(\varepsilon\), already passing matrix construction time for \(\varepsilon = 10^{-2}\).
Regarding overall solve time: already H-LU with the coarsest accuracy needs longer than a simple block-Jacobi preconditioner. Furthermore, matrix-vector solves as are needed for the application of H-LU as a preconditioner are expensive, which already was a problem for Gauss-Seidel preconditioners.
All together, the advantage of H-LU for this specific example will only show for many right-hand sides to solve. As an example, we compare TFQMR+Block-Jacobi (blocksize 100) with Linear Iteration+H-LU ( \(\varepsilon = 10^{-2}\)). Both solvers would need the same overall time for five right-hand side vectors (110.47s vs. 110.73s). Only for more vectors H-LU shows a better performance than Block-Jacobi.
The main problem with block-wise Jacobi for large block sizes was the large setup costs due to inversion of the diagonal matrix blocks. Instead of inversion, these blocks may also be factorised using H-LU factorisation. This permits larger block sizes for block-Jacobi.
𝖧𝖫𝖨𝖡𝗉𝗋𝗈 provides two functions to either restrict a given matrix to it's block-diagonal part
or to copy the block-diagonal of a matrix:
In both cases, the block-diagonal is either defined by the level of the diagonal blocks or by their sizes, e.g., stop recursion if a certain level or a certain blocksize was reached. The default value lvl_leaf signals, that recursion should only be stopped at leaf level (or blocksize = \(n_{\min}\)).
For our problem, a corresponding preconditioner is obtained via
The results for different block sizes are:
blocksize \(100\)
blocksize \(500\)
blocksize \(2000\)
blocksize \(10000\)
blocksize \(20000\)
The results for a blocksize of \(100\) and \(500\) are (almost) identical to the block-Jacobi preconditioner above, only the setup time is reduced, e.g., from 5.19 seconds to 0.58 seconds for blocksize \(500\). The solve times are slightly worse than in the Jacobi version, since matrix-vector solves are involved.
However, already for a blocksize of \(500\), the overall time is faster and even better for a blocksize of \(2000\) with a minimum for \(10000\). Although a larger block size will result in a better convergence but the setup time also increases making block-diagonal H-LU to the fastest solver for our example in terms of overall time.
We have not discussed the accuracy used for H-LU in case of block diagonal preconditioning. For the above tests, the accuracy was \(1\), i.e., a very coarse approximation. Tests with higher accuracy did not resulted in an improvment, for one because the block diagonal structure already sets a lower limit on the approximability of \(A^{-1}\) (see block-Jacobi) and second, because the setup time will then be higher. Only for very large block sizes, this could improve convergence but not necessarily solve times.
Another alternative for computing a H-based preconditioner is to restrict the matrix to the nearfield part only, as defined by the dense matrix blocks. All other low-rank blocks are considered to be zero and will not be updated during LU factorization. This algorithm can be considered as an incomplete LU with limited fill-in.
Since only dense arithmetic is involved, all computations are performed exact w.r.t. machine precision. Furthermore, the computation time is much faster compared to standard H-LU. However, the approximation of \(A^{-1}\) is also limited.
As for block diagonal matrices, we have the following two functions to either restrict a given matrix or to copy the nearfield part:
The second parameter defines, whether low-rank blocks should be deleted/not copied or just zeroed.
In our case, this looks as:
And the results are:
The solve times are not as fast as in the case with block diagonal (Jacobi or H-LU), but not much worse. Especially the fast setup makes this version interesting. Also matrix-vector solves are faster than for standard H-LU.
For the given Helmholtz problem with a single right-hand side, the fastest solver of the block-diagonal H-LU factorisation with a block size of \(10000\). The reason was the fast factorisation of the diagonal blocks, i.e., a low setup time, with a reasonable number of iteration steps.
Using the standard H-LU preconditioner with \(\varepsilon = 10^{-2}\) for comparsion (H-LU time is similar to matrix construction time), it needs 6 right-hand side vectors to be faster than block-diagonal H-LU.
1.8.13