From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 21 Mar 2010 19:22:22 +0000 (+0000)
Subject: working on constraints for matrix blocks
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e22209289c80f569defbc2c54d1dbeea7e277beb;p=dealii-svn.git

working on constraints for matrix blocks

git-svn-id: https://svn.dealii.org/trunk@20871 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/deal.II/include/multigrid/mg_level_object.h b/deal.II/base/include/base/mg_level_object.h
similarity index 97%
rename from deal.II/deal.II/include/multigrid/mg_level_object.h
rename to deal.II/base/include/base/mg_level_object.h
index 1b97b476a8..7b75ed4154 100644
--- a/deal.II/deal.II/include/multigrid/mg_level_object.h
+++ b/deal.II/base/include/base/mg_level_object.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003, 2004, 2005, 2006 by the deal.II authors
+//    Copyright (C) 2003, 2004, 2005, 2006, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -21,9 +21,6 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-/*!@addtogroup mg */
-/*@{*/
-
 /**
  * An array with an object for each level.  The purpose of this class
  * is mostly to store objects and allow access by level number, even
@@ -34,7 +31,9 @@ DEAL_II_NAMESPACE_OPEN
  * In most cases, the objects which are stored on each levels, are
  * either matrices or vectors.
  *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2005
+ * @ingroup mg
+ * @ingroup data
+ * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2005, 2010
  */
 template<class Object>
 class MGLevelObject : public Subscriptor
@@ -116,7 +115,6 @@ class MGLevelObject : public Subscriptor
     std::vector<std_cxx1x::shared_ptr<Object> > objects;
 };
 
-/*@}*/
 
 /* ------------------------------------------------------------------- */
 
diff --git a/deal.II/deal.II/include/multigrid/mg_block_smoother.h b/deal.II/deal.II/include/multigrid/mg_block_smoother.h
index 3065fc7e90..b51afecf49 100644
--- a/deal.II/deal.II/include/multigrid/mg_block_smoother.h
+++ b/deal.II/deal.II/include/multigrid/mg_block_smoother.h
@@ -20,7 +20,7 @@
 #include <lac/vector_memory.h>
 #include <lac/block_vector.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
diff --git a/deal.II/deal.II/include/multigrid/mg_matrix.h b/deal.II/deal.II/include/multigrid/mg_matrix.h
index d66cf2c3de..274b1eb3ea 100644
--- a/deal.II/deal.II/include/multigrid/mg_matrix.h
+++ b/deal.II/deal.II/include/multigrid/mg_matrix.h
@@ -16,7 +16,7 @@
 #include <lac/vector.h>
 #include <lac/sparse_matrix.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 DEAL_II_NAMESPACE_OPEN
 
diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h
index d298470a95..e28f7cd486 100644
--- a/deal.II/deal.II/include/multigrid/mg_smoother.h
+++ b/deal.II/deal.II/include/multigrid/mg_smoother.h
@@ -19,7 +19,7 @@
 #include <lac/pointer_matrix.h>
 #include <lac/vector_memory.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <vector>
 
 DEAL_II_NAMESPACE_OPEN
diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.h b/deal.II/deal.II/include/multigrid/mg_transfer.h
index 8c63dc237a..7db0779cc5 100644
--- a/deal.II/deal.II/include/multigrid/mg_transfer.h
+++ b/deal.II/deal.II/include/multigrid/mg_transfer.h
@@ -27,7 +27,7 @@
 #include <lac/vector_memory.h>
 
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 
diff --git a/deal.II/deal.II/include/multigrid/mg_transfer_block.h b/deal.II/deal.II/include/multigrid/mg_transfer_block.h
index bc895bc3ec..a9f7451ef2 100644
--- a/deal.II/deal.II/include/multigrid/mg_transfer_block.h
+++ b/deal.II/deal.II/include/multigrid/mg_transfer_block.h
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors
+//    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -27,7 +27,7 @@
 #include <lac/block_matrix_array.h>
 
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 
diff --git a/deal.II/deal.II/include/multigrid/mg_transfer_component.h b/deal.II/deal.II/include/multigrid/mg_transfer_component.h
index 2d1e1de352..3a209635f4 100644
--- a/deal.II/deal.II/include/multigrid/mg_transfer_component.h
+++ b/deal.II/deal.II/include/multigrid/mg_transfer_component.h
@@ -27,7 +27,7 @@
 #include <lac/vector_memory.h>
 
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 
diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h
index 8442d467ef..fbcccdce6b 100644
--- a/deal.II/deal.II/include/multigrid/multigrid.h
+++ b/deal.II/deal.II/include/multigrid/multigrid.h
@@ -21,7 +21,7 @@
 #include <lac/sparse_matrix.h>
 #include <lac/vector.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <multigrid/mg_dof_handler.h>
 
 #include <vector>
diff --git a/deal.II/deal.II/include/numerics/mesh_worker.h b/deal.II/deal.II/include/numerics/mesh_worker.h
index 3602b4331f..6f9eda3ae4 100644
--- a/deal.II/deal.II/include/numerics/mesh_worker.h
+++ b/deal.II/deal.II/include/numerics/mesh_worker.h
@@ -18,7 +18,7 @@
 #include <base/geometry_info.h>
 #include <base/named_data.h>
 #include <lac/block_indices.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <numerics/mesh_worker_vector_selector.h>
 #include <numerics/mesh_worker_info.h>
 #include <numerics/mesh_worker_assembler.h>
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h
index 01371224c7..e1af235b0b 100644
--- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h
+++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h
@@ -17,7 +17,7 @@
 #include <base/named_data.h>
 #include <base/smartpointer.h>
 #include <lac/block_vector.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 DEAL_II_NAMESPACE_OPEN
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_output.h b/deal.II/deal.II/include/numerics/mesh_worker_output.h
index 0e557c67af..1ae30f1ec8 100644
--- a/deal.II/deal.II/include/numerics/mesh_worker_output.h
+++ b/deal.II/deal.II/include/numerics/mesh_worker_output.h
@@ -17,7 +17,7 @@
 #include <base/named_data.h>
 #include <base/smartpointer.h>
 #include <lac/block_vector.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 
 DEAL_II_NAMESPACE_OPEN
diff --git a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
index 80312e8f88..41f709fd65 100644
--- a/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
+++ b/deal.II/deal.II/include/numerics/mesh_worker_vector_selector.h
@@ -16,7 +16,7 @@
 #include <base/named_data.h>
 #include <base/tensor.h>
 #include <base/smartpointer.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 DEAL_II_NAMESPACE_OPEN
 
diff --git a/deal.II/deal.II/source/multigrid/mg_base.cc b/deal.II/deal.II/source/multigrid/mg_base.cc
index 90476e759b..e7088cdfd7 100644
--- a/deal.II/deal.II/source/multigrid/mg_base.cc
+++ b/deal.II/deal.II/source/multigrid/mg_base.cc
@@ -17,7 +17,7 @@
 #include <fe/fe.h>
 #include <multigrid/mg_dof_handler.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <multigrid/mg_tools.h>
 
 
diff --git a/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc b/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc
index e85e42f8fd..8a7fbba6f1 100644
--- a/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc
+++ b/deal.II/deal.II/source/multigrid/mg_tools.all_dimensions.cc
@@ -18,7 +18,7 @@
 #include <lac/block_vector.h>
 #include <multigrid/mg_tools.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 
 #include <vector>
 #include <algorithm>
diff --git a/deal.II/deal.II/source/multigrid/mg_tools.cc b/deal.II/deal.II/source/multigrid/mg_tools.cc
index 92eaf98f41..035e80a42c 100644
--- a/deal.II/deal.II/source/multigrid/mg_tools.cc
+++ b/deal.II/deal.II/source/multigrid/mg_tools.cc
@@ -26,7 +26,7 @@
 #include <multigrid/mg_dof_accessor.h>
 #include <multigrid/mg_tools.h>
 #include <multigrid/mg_base.h>
-#include <multigrid/mg_level_object.h>
+#include <base/mg_level_object.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe.h>
 
diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc
index 8eff89bddf..64322d2285 100644
--- a/deal.II/examples/step-12/step-12.cc
+++ b/deal.II/examples/step-12/step-12.cc
@@ -456,14 +456,14 @@ void Step12<dim>::integrate_cell_term (DoFInfo& dinfo, CellInfo& info)
 }
 
 				 // Now the same for the boundary terms. Note
-				 // that now we use FEFaceValuesBase, the base
+				 // that now we use FEValuesBase, the base
 				 // class for both FEFaceValues and
 				 // FESubfaceValues, in order to get access to
 				 // normal vectors.
 template <int dim>
 void Step12<dim>::integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info)
 {
-  const FEFaceValuesBase<dim>& fe_v = info.fe_values();
+  const FEValuesBase<dim>& fe_v = info.fe_values();
   FullMatrix<double>& local_matrix = dinfo.matrix(0).matrix;
   Vector<double>& local_vector = dinfo.vector(0).block(0);
 
@@ -512,14 +512,14 @@ void Step12<dim>::integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2,
 {
 				   // For quadrature points, weights,
 				   // etc., we use the
-				   // FEFaceValuesBase object of the
+				   // FEValuesBase object of the
 				   // first argument.
-  const FEFaceValuesBase<dim>& fe_v = info1.fe_values();
+  const FEValuesBase<dim>& fe_v = info1.fe_values();
 
 				   // For additional shape functions,
 				   // we have to ask the neighbors
-				   // FEFaceValuesBase.
-  const FEFaceValuesBase<dim>& fe_v_neighbor = info2.fe_values();
+				   // FEValuesBase.
+  const FEValuesBase<dim>& fe_v_neighbor = info2.fe_values();
 
 				   // Then we get references to the
 				   // four local matrices. The letters
diff --git a/deal.II/lac/include/lac/constraint_matrix.h b/deal.II/lac/include/lac/constraint_matrix.h
index c4c21c4727..79ffcfe2b3 100644
--- a/deal.II/lac/include/lac/constraint_matrix.h
+++ b/deal.II/lac/include/lac/constraint_matrix.h
@@ -1825,6 +1825,18 @@ class ConstraintMatrix : public Subscriptor
 				 const Table<2,bool>             &dof_mask,
 				 internal::bool2type<true>) const;
 
+				     /**
+				      * Internal helper function for
+				      * distribute_local_to_global
+				      * function.
+				      *
+				      * Creates a list of affected
+				      * rows for sistribution.
+				      */
+    void
+    make_sorted_row_list (const std::vector<unsigned int> &local_dof_indices,
+			  internals::GlobalRowsFromLocal  &global_rows) const;
+
 				     /**
 				      * Internal helper function for
 				      * distribute_local_to_global function.
diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h
index 5f20f9465f..d9a2516ad3 100644
--- a/deal.II/lac/include/lac/constraint_matrix.templates.h
+++ b/deal.II/lac/include/lac/constraint_matrix.templates.h
@@ -1722,6 +1722,85 @@ namespace internals
 } // end of namespace internals
 
 
+void
+ConstraintMatrix::
+make_sorted_row_list (const std::vector<unsigned int> &local_dof_indices,
+		      internals::GlobalRowsFromLocal  &global_rows) const
+{
+  const unsigned int n_local_dofs = local_dof_indices.size();
+				   // when distributing the local data to
+				   // the global matrix, we can quite
+				   // cheaply sort the indices (obviously,
+				   // this introduces the need for
+				   // allocating some memory on the way, but
+				   // we need to do this only for rows,
+				   // whereas the distribution process
+				   // itself goes over rows and
+				   // columns). This has the advantage that
+				   // when writing into the global matrix,
+				   // we can make use of the sortedness.
+
+				   // so the first step is to create a
+				   // sorted list of all row values that are
+				   // possible. these values are either the
+				   // rows from unconstrained dofs, or some
+				   // indices introduced by dofs constrained
+				   // to a combination of some other
+				   // dofs. regarding the data type, choose
+				   // an STL vector of a pair of unsigned
+				   // ints (for global columns) and internal
+				   // data (containing local columns +
+				   // possible jumps from
+				   // constraints). Choosing an STL map or
+				   // anything else M.K. knows of would be
+				   // much more expensive here!
+
+				   // cache whether we have to resolve any
+				   // indirect rows generated from resolving
+				   // constrained dofs.
+  unsigned int added_rows = 0;
+
+				   // first add the indices in an unsorted
+				   // way and only keep track of the
+				   // constraints that appear. They are
+				   // resolved in a second step.
+  for (unsigned int i = 0; i<n_local_dofs; ++i)
+    {
+      if (is_constrained(local_dof_indices[i]) == false)
+	{
+	  global_rows.global_row(added_rows)  = local_dof_indices[i];
+	  global_rows.local_row(added_rows++) = i;
+	  continue;
+	}
+      global_rows.insert_constraint(i);
+    }
+  global_rows.sort();
+
+  const unsigned int n_constrained_rows = n_local_dofs-added_rows;
+  for (unsigned int i=0; i<n_constrained_rows; ++i)
+    {
+      const unsigned int local_row = global_rows.last_constrained_local_row();
+      Assert (local_row < n_local_dofs, ExcIndexRange(local_row, 0, n_local_dofs));
+      const unsigned int global_row = local_dof_indices[local_row];
+      Assert (is_constrained(global_row), ExcInternalError());
+      const ConstraintLine & position =
+	lines[lines_cache[calculate_line_index(global_row)]];
+      if (position.inhomogeneity != 0)
+	global_rows.set_last_row_inhomogeneous();
+      else
+	global_rows.pop_back();
+      for (unsigned int q=0; q<position.entries.size(); ++q)
+	global_rows.insert_index (position.entries[q].first,
+				  local_row,
+				  position.entries[q].second);
+    }
+}
+
+
+
+
+//TODO: This function is DANGEROUS, because it does more than it claims!
+
 
 				// Basic idea of setting up a list of
 				// all global dofs: first find all rows and columns
diff --git a/deal.II/lac/include/lac/matrix_block.h b/deal.II/lac/include/lac/matrix_block.h
index 7d251af375..668fa620cd 100644
--- a/deal.II/lac/include/lac/matrix_block.h
+++ b/deal.II/lac/include/lac/matrix_block.h
@@ -27,30 +27,52 @@ DEAL_II_NAMESPACE_OPEN
  * A wrapper around a matrix object, storing the coordinates in a
  * block matrix as well.
  *
- * @note They are the position in the global matrix, where we don't use
- * BlockSparseMatrix. Matrices may have the same coordinates for the
- * following reason: A preconditioner for the Oseen system can be
- * built as a block system, where the pressure block is of the form
- * M^-1 F A^-1 with m the mass matrix (u,v), A the Laplacian and F the
- * advection diffusion operator. So, I need to build three matrices
- * for a single block in my system + in some cases a pressure
- * stabilization.
+ * This class is an alternative to BlockMatrixBase, if you only want
+ * to generate a single block of the system, not the whole
+ * system. Using the add() functions of this class, it is possible to
+ * use the standard assembling functions used for block matrices, but
+ * only enter in one of the blocks and still avoiding the index
+ * computations innvolved.
+
+ * The reason for this class is, that we may need a different number
+ * of matrices for different blocks in a block system. For example, a
+ * preconditioner for the Oseen system can be built as a block system,
+ * where the pressure block is of the form
+ * <b>M</b><sup>-1</sup><b>FA</b><sup>-1</sup> with <b>M</b> the
+ * pressure mass matrix, <b>A</b> the pressure Laplacian and <b>F</b>
+ * the advection diffusion operator applied to the pressure
+ * space. Since only a single matrix is needed for the other blocks,
+ * using BlockSparseMatrix or similar would be a waste of memory.
+ *
+ * While the add() functions make a MatrixBlock appear like a block
+ * matrix for assembling, the functions vmult(), Tvmult(),
+ * vmult_add(), and Tvmult_add() make it behave like a MATRIX, when it
+ * comes to applying it to a vector. This behavior allows us to store
+ * MatrixBlock objects in vectors, for instance in MGLevelObject
+ * without extracting the #matrix first. 
  *
  * MatrixBlock comes handy when using BlockMatrixArray. Once the
  * MatrixBlock has been properly initalized and filled, it can be used
- * as:
+ * in the simplest case as:
  * @code
- * std::vector<MatrixBlock<SparseMatrix<double> > > blocks;
+ * MatrixBlockVector<SparseMatrix<double> > > blocks;
  *
  * ...
  *
  * BlockMatrixArray matrix (n_blocks, n_blocks);
  *
  * for (unsigned int i=0;i<blocks.size;++i)
- *   matrix.enter(blocks[i].row, blocks[i].column, blocks[i].matrix);
+ *   matrix.enter(blocks.block(i).row, blocks.block(i).column, blocks.matrix(i));
  * @endcode
  *
+ * Here, we have not gained very much, except that we do not need to
+ * set up empty blocks in the block system.
+ *
+ * @todo Example for the product preconditioner of the pressure Schur
+ * complement.
  *
+ * @ingroup Matrix2
+ * @ingroup vector_valued
  * @author Guido Kanschat, 2006
  */
 template <class MATRIX>
@@ -79,33 +101,34 @@ class MatrixBlock
 
 				     /**
 				      * Add <tt>value</tt> to the
-				      * element (<i>i,j</i>).  Throws
+				      * element (<i>i,j</i>). Throws
 				      * an error if the entry does not
-				      * exist or if <tt>value</tt> is
-				      * not a finite number. Still, it
-				      * is allowed to store zero
-				      * values in non-existent fields.
+				      * exist or if it is in a
+				      * different block.
 				      */
     void add (const unsigned int i,
               const unsigned int j,
 	      const typename MATRIX::value_type value);
 
                                        /**
-                                        * Add all elements given in a
-                                        * FullMatrix<double> into sparse
+                                        * Add all elements in a
+                                        * FullMatrix into sparse
                                         * matrix locations given by
-                                        * <tt>indices</tt>. In other words,
-                                        * this function adds the elements in
-                                        * <tt>full_matrix</tt> to the
-                                        * respective entries in calling
-                                        * matrix, using the local-to-global
-                                        * indexing specified by
-                                        * <tt>indices</tt> for both the rows
-                                        * and the columns of the
-                                        * matrix. This function assumes a
-                                        * quadratic sparse matrix and a
-                                        * quadratic full_matrix, the usual
-                                        * situation in FE calculations.
+                                        * <tt>indices</tt>. This function
+                                        * assumes a quadratic sparse
+                                        * matrix and a quadratic
+                                        * full_matrix.  The global
+                                        * locations are translated
+                                        * into locations in this block
+                                        * and ExcBlockIndexMismatch is
+                                        * thrown, if the global index
+                                        * does not point into the
+                                        * block refered to by #row and
+                                        * #column.
+					*
+					* @todo
+					* <tt>elide_zero_values<tt> is
+					* currently ignored.
 					*
 					* The optional parameter
 					* <tt>elide_zero_values</tt> can be
@@ -123,11 +146,33 @@ class MatrixBlock
 	      const bool                       elide_zero_values = true);
 
                                        /**
-                                        * Same function as before, but now
-                                        * including the possibility to use
-                                        * rectangular full_matrices and
-                                        * different local-to-global indexing
-                                        * on rows and columns, respectively.
+					* Add all elements in a
+                                        * FullMatrix into global
+                                        * locations given by
+                                        * <tt>row_indices</tt> and
+                                        * <tt>col_indices</tt>,
+                                        * respectively. The global
+                                        * locations are translated
+                                        * into locations in this block
+                                        * and ExcBlockIndexMismatch is
+                                        * thrown, if the global index
+                                        * does not point into the
+                                        * block refered to by #row and
+                                        * #column.
+					*
+					* @todo
+					* <tt>elide_zero_values<tt> is
+					* currently ignored.
+					*
+					* The optional parameter
+					* <tt>elide_zero_values</tt> can be
+					* used to specify whether zero
+					* values should be added anyway or
+					* these should be filtered away and
+					* only non-zero data is added. The
+					* default value is <tt>true</tt>,
+					* i.e., zero values won't be added
+					* into the matrix.
 					*/
     template <typename number>
     void add (const std::vector<unsigned int>& row_indices,
@@ -137,10 +182,27 @@ class MatrixBlock
 
                                        /**
                                         * Set several elements in the
-                                        * specified row of the matrix with
-                                        * column indices as given by
-                                        * <tt>col_indices</tt> to the
-                                        * respective value.
+                                        * specified row of the matrix
+                                        * with column indices as given
+                                        * by <tt>col_indices</tt> to
+                                        * the respective value. This
+                                        * is the function doing thye
+                                        * actual work for the ones
+                                        * adding full matrices. The
+                                        * global locations
+                                        * <tt>row_index</tt> and
+                                        * <tt>col_indices</tt> are
+                                        * translated into locations in
+                                        * this block and
+                                        * ExcBlockIndexMismatch is
+                                        * thrown, if the global index
+                                        * does not point into the
+                                        * block refered to by #row and
+                                        * #column.
+					*
+					* @todo
+					* <tt>elide_zero_values<tt> is
+					* currently ignored.
 					*
 					* The optional parameter
 					* <tt>elide_zero_values</tt> can be
@@ -153,7 +215,7 @@ class MatrixBlock
 					* into the matrix.
 					*/
     template <typename number>
-    void add (const unsigned int               row,
+    void add (const unsigned int               row_index,
 	      const std::vector<unsigned int>& col_indices,
 	      const std::vector<number>&       values,
 	      const bool                       elide_zero_values = true);
@@ -183,6 +245,50 @@ class MatrixBlock
 	      const bool          elide_zero_values = true,
 	      const bool          col_indices_are_sorted = false);
 
+				     /**
+				      * Matrix-vector-multiplication,
+				      * forwarding to the same
+				      * function in MATRIX. No index
+				      * computations are done, thus,
+				      * the vectors need to have sizes
+				      * matching #matrix.
+				      */
+    template<class VECTOR>
+    void vmult (VECTOR& w, const VECTOR& v) const;
+    
+				     /**
+				      * Matrix-vector-multiplication,
+				      * forwarding to the same
+				      * function in MATRIX. No index
+				      * computations are done, thus,
+				      * the vectors need to have sizes
+				      * matching #matrix.
+				      */
+    template<class VECTOR>
+    void vmult_add (VECTOR& w, const VECTOR& v) const;
+    
+				     /**
+				      * Matrix-vector-multiplication,
+				      * forwarding to the same
+				      * function in MATRIX. No index
+				      * computations are done, thus,
+				      * the vectors need to have sizes
+				      * matching #matrix.
+				      */
+    template<class VECTOR>
+    void Tvmult (VECTOR& w, const VECTOR& v) const;
+
+				     /**
+				      * Matrix-vector-multiplication,
+				      * forwarding to the same
+				      * function in MATRIX. No index
+				      * computations are done, thus,
+				      * the vectors need to have sizes
+				      * matching #matrix.
+				      */
+    template<class VECTOR>
+    void Tvmult_add (VECTOR& w, const VECTOR& v) const;
+
 				     /**
 				      * The block number computed from
 				      * an index by using
@@ -229,6 +335,8 @@ class MatrixBlock
  * pointers, in order to allow for copying and rearranging. Each
  * matrix block can be identified by name.
  *
+ * @relates MatrixBlock
+ * @ingroup vector_valued
  * @author Baerbel Janssen, Guido Kanschat, 2010
  */
 template <class MATRIX>
@@ -264,6 +372,48 @@ class MatrixBlockVector : public NamedData<boost::shared_ptr<MatrixBlock<MATRIX>
 };
 
 
+/**
+ * A vector of MGLevelObject<MatrixBlock>, which is implemented using shared
+ * pointers, in order to allow for copying and rearranging. Each
+ * matrix block can be identified by name.
+ *
+ * @relates MatrixBlock
+ * @ingroup vector_valued
+ * @author Baerbel Janssen, Guido Kanschat, 2010
+ */
+template <class MATRIX>
+class MGMatrixBlockVector : public NamedData<boost::shared_ptr<MatrixBlock<MATRIX> > >
+{
+  public:
+				     /**
+				      * Add a new matrix block at the
+				      * position <tt(row,column)</tt>
+				      * in the block system.
+				      */
+    void add(unsigned int row, unsigned int column,
+	     const std::string& name,
+	     const BlockIndices* block_indices);
+				     /**
+				      * Access a constant reference to
+				      * the block at position <i>i</i>.
+				      */
+    const MatrixBlock<MATRIX>& block(unsigned int i) const;
+				     /**
+				      * Access a reference to
+				      * the block at position <i>i</i>.
+				      */
+    MatrixBlock<MATRIX>& block(unsigned int i);
+				     /**
+				      * Access the matrix at position
+				      * <i>i</i> for read and write
+				      * access.
+				      */
+    MATRIX& matrix(unsigned int i);
+    
+    
+};
+
+
 //----------------------------------------------------------------------//
 
 template <class MATRIX>
@@ -336,6 +486,43 @@ MatrixBlock<MATRIX>::add (const std::vector<unsigned int>&         row_indices,
 }
 
 
+template <class MATRIX>
+template <typename number>
+inline
+void
+MatrixBlock<MATRIX>::add (const unsigned int   b_row,
+			  const unsigned int   n_cols,
+			  const unsigned int  *col_indices,
+			  const number        *values,
+			  const bool,
+			  const bool)
+{
+  Assert(block_indices != 0, ExcNotInitialized());
+  const std::pair<unsigned int, unsigned int> bi
+    = block_indices->global_to_local(b_row);
+
+				   // In debug mode, we check whether
+				   // all indices are in the correct
+				   // block.
+
+				   // Actually, for the time being, we
+				   // leave it at this. While it may
+				   // not be the most efficient way,
+				   // it is at least thread safe.
+//#ifdef DEBUG
+  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
+
+  for (unsigned int j=0;j<n_cols;++j)
+    {
+      const std::pair<unsigned int, unsigned int> bj
+	= block_indices->global_to_local(col_indices[j]);
+      Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
+      
+      matrix.add(bi.second, bj.second, values[j]);
+    }
+//#endif
+}
+
 
 template <class MATRIX>
 template <typename number>
@@ -374,45 +561,46 @@ MatrixBlock<MATRIX>::add (const unsigned int               row,
 
 
 template <class MATRIX>
-template <typename number>
+template <class VECTOR>
 inline
 void
-MatrixBlock<MATRIX>::add (const unsigned int   b_row,
-			  const unsigned int   n_cols,
-			  const unsigned int  *col_indices,
-			  const number        *values,
-			  const bool,
-			  const bool)
+MatrixBlock<MATRIX>::vmult (VECTOR& w, const VECTOR& v) const
 {
-  Assert(block_indices != 0, ExcNotInitialized());
-  const std::pair<unsigned int, unsigned int> bi
-    = block_indices->global_to_local(b_row);
+  matrix.vmult(w,v);
+}
 
-				   // In debug mode, we check whether
-				   // all indices are in the correct
-				   // block.
 
-				   // Actually, for the time being, we
-				   // leave it at this. While it may
-				   // not be the most efficient way,
-				   // it is at least thread safe.
-//#ifdef DEBUG
-  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
+template <class MATRIX>
+template <class VECTOR>
+inline
+void
+MatrixBlock<MATRIX>::vmult_add (VECTOR& w, const VECTOR& v) const
+{
+  matrix.vmult_add(w,v);
+}
 
-  for (unsigned int j=0;j<n_cols;++j)
-    {
-      const std::pair<unsigned int, unsigned int> bj
-	= block_indices->global_to_local(col_indices[j]);
-      Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
-      
-      matrix.add(bi.second, bj.second, values[j]);
-    }
-  
-//#endif
 
+template <class MATRIX>
+template <class VECTOR>
+inline
+void
+MatrixBlock<MATRIX>::Tvmult (VECTOR& w, const VECTOR& v) const
+{
+  matrix.Tvmult(w,v);
+}
+
+
+template <class MATRIX>
+template <class VECTOR>
+inline
+void
+MatrixBlock<MATRIX>::Tvmult_add (VECTOR& w, const VECTOR& v) const
+{
+  matrix.Tvmult_add(w,v);
 }
 
 
+
 //----------------------------------------------------------------------//
 
 template <class MATRIX>