- Generated by
1.9.8
|
HLIBpro 3.1
|
For equation systems with many right hand sides, e.g.
\[ A x_i = b_i \]
with \(i = 1 \ldots n_{rhs}, n_{rhs} \gg 1,\) different ways are supported in π§π«π¨π‘πππ to solve the equations, which shall be introduced in this example.
We assume a given system matrix \(A\) and right hand sides (RHSs) \(b_i\). We are looking for the solutions \(x_i\) in \( A x_i = b_i \).
As with all previous applications of H-matrices, the H-LU factorisation shall be used to solve the equation systems. This may either be performed directly, i.e., as a direct solver involving a single matrix-vector evaluation, or with the help of an iterative solver. In all examples so far, the iterative approach was used since it usually relaxes the accuracy demands on the H-LU factorisation and hence reduces the factorisation runtime.
However, iterative solvers perform several matrix-vector multiplications for a single RHS. Therefore, this approach may increase the overall time significantly for many RHSs.
Alternatively, the direct solver approach may be used, which normally demands a much higher accuracy for H-LU to obtain the corresponding accuracy in the solution vector. A higher H-LU accuracy often increases the runtime and memory requirement of the factorisation process. However, depending on the number of RHSs, this additional costs may shorten the overall solution process.
Beside solving for a single RHS, π§π«π¨π‘πππ also offers to solve for all RHSs simultaneously when used as a direct solver, as described below.
In the following, we use a system matrix \(A\) and corresponding RHSs \(b_i\) from the Boundary Element Methods example. However, instead of the Helmholtz kernel, the Laplace kernel is used. Furthermore, the matrix is constructed as an unsymmetric matrix, i.e., with upper and lower triangular part (see below for standard symmetric case).
The initialisation and construction of cluster trees is:
This is followed by H-matrix construction:
Here, the unsymmetric matrix is constructed instead of the standard symmetric version.
Next, the RHSs shall be computed. For the BEM problem this is accomplished by deriving from TBEMFunction. Since multiple, different RHSs are to be constructed, two parameters alpha and beta are added, which are used in the evaluation function:
Up to now, the RHS was represented as a single vector. For multiple RHSs this may be extended to a list of vectors. However, even more compact, and also useful for the different solve methods, is the representation as a dense matrix, i.e., each column corresponds to a different RHS:
\[ B := \left[ b_0 b_1 \ldots \right] \]
In the source code the matrix \(B \in R^{n \times n_{rhs}}\) is stored in the dense matrix RHS:
For the actual construction, each column of RHS is taken as a standard vector and computed in parallel using tbb::parallel_for :
The parameters alpha and beta of TMyRHS are chosen randomly. Similar to vectors, also the dense matrix has to be reordered to have the same ordering as the indices in the H-matrix. The permutation is here only applied to the rows, since they correspond to geometrical indices, while the column index only represents the index of the RHS.
TScalarVector, which should be true for almost all cases.For all methods to solve the above equations, the LU factorisation \(A = LU\) is needed, which is computed next:
First, the classical method using an iterative solver is used. For this, each equation \(A x_i = b_i\) is solved in parallel using the solve function. Analogously to the RHSs, the solution vectors are stored in a dense matrix SOL:
After all equations are solved, the residual is computed with respect to all RHSs:
Next, the \H-LU should be used as a direct solver. The only difference to the previous method is, that instead of using the function solve, the function apply of the operator A_inv is used, which applies the inverse of \(A\) to the RHS \(b_i\):
Compared to the iterative method, this approach is faster but also the residual is usually higher. To achieve the same approximation as with the previous method, the H-LU factorisation has to be performed with an higher accuracy.
Finally, the LU factors \(L\) and \(U\) are used directly to solve all equations at once. Due to
\[A X = L U X = B\]
with \(X \in R^{n \times n_{rhs}}\) being the solution matrix, first with the lower triangular matrix \(L\) is solved and afterwards with the upper triangular matrix \(U\):
Since \(L\) is computed with unit diagonal, this also has to be defined for solving. Furthermore, instead of A_inv, which also represents the inverse operator, the matrix A_copy has to be used since it holds the actual LU factors.
apply_adjoint is to be used for the second solve_lower_left call.The matrix approach yields to same residual as in the previous case. However, it may be more efficient in terms of runtime than the vector approach since it can exploit data locality much better.
However, the direct vector approach scales better with the number of cores. Hence, depending on your CPU configuration, this approach may be faster. Of course, the number of RHSs should be reasonably large for a good parallel scaling.
1.9.8