Fix: Cope with omitted diagonal blocks in outer GS loop.
The former implementation had two shortcomings:
- It assumed the diagonal block to exist. This is not necessarily the case, e.g., in the case of a
BCRSMatrixwhere a non-retrievable block within the dimensions of the matrix is equivalent to it being zero. Nevertheless we should be able to cope with this case, in particular because semi-definite matrices should be covered by this approach due to the inherent regularization. - Access to the diagonal element was performed via the random access operator which is relatively costly.
The new implementation tackles both problems. In case the diagonal block does not exist, a temporary block is constructed and passed to the local solver instead.
On a side note: A more generic approach to creating and initializing the temporary block might be needed for more complex nested matrix schemes.
Merge request reports
Activity
Added 1 commit:
- 4ff9dbfb - Fix: Cope with omitted diagonal blocks in outer GS loop.
Added 1 commit:
- bfeb2942 - Fix: Cope with omitted diagonal blocks in outer GS loop.
46 using Block = typename M::block_type; 47 const Block* diag = nullptr; 48 for (auto cIt = row_i.begin(); cIt != row_i.end(); ++cIt) { 49 size_t j = cIt.index(); 50 cIt->mmv(x[j], ri); 51 if (j == i) 52 diag = &*cIt; 53 } 47 54 48 55 std::bitset<V::block_type::dimension> ignore_i(0); 49 56 if (ignore != nullptr) 50 57 ignore_i = (*ignore)[i]; 58 51 59 // Update iterate with correction 52 x[i] += localSolver(row_i[i], std::move(ri), ignore_i); 60 x[i] += localSolver(diag ? *diag : Block(0.0), std::move(ri), ignore_i); Added 2 commits:
Mentioned in commit c1694ae1
What would you propose to do instead?
I think that for now it is at least consistent in the sense that a non-existing block is treated the same way as a zero-block (because that's what a non-existing block is supposed to represent anyway imho).
The cases that come to my mind are the following:
- The ignore nodes are set for the entire block: Then the block correction should be zero in any case (and the available ignore node handling does so). I though of the small instantiation overhead to be negligible.
- The ignore nodes are not set for the entire block: Then something should be computed, although I'm honestly not too confident about the right way to handle this generically. Until now I thought of this problem as being use-case specific and whoever is applying a BlockGS to such a matrix should be aware of what he expects to happen and eventually modify his local solving routines accordingly. In this case however, it affects the outer GS loop ...
I am very willing to do further tweaking here, but I also want to keep clarity about what is happening conceptually and what are the distinct roles of ignore nodes and indefinite blocks.
I'd say that conceptually, the local solver should compute an approximation for the local problem. Using approximations instead of exact solutions is only OK since we compute corrections. Since we also want to deal with spsd matrices we allow local solvers aware of this situation. The aim of the local solver is then to compute an approximation of the application of the pseudo-monotone of the diagonal block. In case of a zero diagonal block the pseudo monotone is also zero an the only reasonable approximation to the solution would be zero. For anything else it would be very hard to show convergence.
-
From the discussion we just had, I understood that the current implementation does the right thing in the sense that for zero blocks the residual is assumed to be appropriate anyways (i.e. has to be zero). However
- there is overhead in instantiating a temporary for something that is expected to compute a zero correction and
- the current way of explicitly instantiating a zero block is not suited for generic matrix-nesting schemes.
I added related todos to #8 and will tackle them in the future.
