<div class="gmail_quote">On Fri, May 20, 2011 at 15:16, Thomas Witkowski <span dir="ltr"><<a href="mailto:thomas.witkowski@tu-dresden.de">thomas.witkowski@tu-dresden.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">
Thanks for the information. One related question: The system that is solved here is a 3x3 block matrix: Two of these blocks are empty, two are the discretized (FEM) identity, three are the discretized Laplace operator</blockquote>
<div><br></div><div>What equation? Have you assumed special boundary conditions? Problems like elasticity and Stokes normally need a symmetric stress tensor. Only for constant coefficients and certain boundary conditions, is it okay to use the "block Laplacian" form.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"> and one block is the discretization of a complex sum of operators. At all, there are only three different blocks. Can I make somehow advantages of this structure? I could create the three inverses of the three different blocks, but exists than a cheap/simple way to produce the action of the inverse of the whole 3x3 block matrix?</blockquote>
</div><br><div>For the special case of a Kronecker product, you can solve using solves with the blocks, but not in the general case. Consider the structured grid 2D Laplacian [A, -I; -I, A] where A is the tridiagonal problem for a 1D Laplacian plus 2*I on the diagonal. If you could build a direct solver for the 2D using direct solvers for the 1D problem, then you would have an O(n) algorithm for a direct solve of a 2D Laplacian.</div>
<div><br></div><div>You can, however, build an iterative method from the blocks. If each block is built separately, you can (with petsc-dev) put them into MatNest and use PCFieldSplit with direct solves in splits.</div><div>
<br></div><div>An alternative is to interlace the degrees of freedom and use the (S)BAIJ matrix format. This will reduce the amount of memory needed for your matrix (by using a more efficient storage format) and will improve the memory performance. Unfortunately, there are no third-party direct solvers for BAIJ. Note that your assembly code can always call MatSetValuesBlocked() (or ..BlockedLocal()) regardless of the matrix format, so you can switch between AIJ and BAIJ at run-time.</div>
<div><br></div>Why bother with a direct solver if your problem has constant coefficients?