Glas List: [glas] (no subject)

From: Wolfgang Bangerth (bangerth_at_[hidden])
Date: 2005-10-13 19:01:50

In reply to an already several days old message:

> I recall a demand for supporting vectors of matrices and matrices of
> matrices. My opinion is that this should be another type of collection
> where the value_type is still the value_type of the lowest level of
> matrix in the nesting of matrices.

There is definitely a demand for something like this, though I wouldn't 
suggest that implementing them as matrices of matrices is a particularly 
clever way to go.

The application where this is important is when one can split a matrix 
into blocks. Typical applications are discretizations of vector-valued 
partial differential equations (think: mixed methods, saddle point 
problems) where instead of solving with a matrix directly, one would like 
to form Schur complements of its blocks. The right abstraction therefore 
is to consider the whole matrix as an array of submatrices (blocks).

The way we implement this in deal.II (and I believe this is the correct 
way) is indeed as an array of matrices. One can access individual blocks, 
but if you try to access individual elements of the matrix, one gets a 
reference to the corresponding element in the block in which it resides. 
The BlockMatrix class therefore has a dual interface: access the blocks of 
the array through
   Matrix & block (const unsigned int block_i, const unsigned int block_j)
and access the elements of the matrix through something like
   reference operator() (const unsigned int i, const unsigned int j) {
     return blocks[block_from_global(i)][block_from_global(j)]
             ->operator() (local_from_global(i), local_from_global(j));
   }

Through the second interface, one can hide the actual storage scheme from 
algorithms that only expect a matrix, while the first interface allows to 
extract and operate on individual blocks.

If, on the contrary, this would be implemented as matrices-of-matrices, 
then one would only have the first interface, but algorithms that just 
expect a transparent matrix object won't work because your operator() 
doesn't return individual entries of the matrix, but an entire block.

For an example of an implementation of such concepts, feel free to take a 
look at
   http://www.dealii.org/5.2.0/doxygen/lac/classBlockSparseMatrix.html

Best
   Wolfgang

-------------------------------------------------------------------------
Wolfgang Bangerth                email:            bangerth_at_[hidden]
                                  www: http://www.math.tamu.edu/~bangerth/