From 26c34f6b45bdd540aaf70c158ca3f4653d3e8a7d Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 22 Jun 2001 10:36:46 +0000
Subject: [PATCH] New class CompressedSparsityPattern. Use it. New test.

git-svn-id: https://svn.dealii.org/trunk@4749 0785d39b-7218-0410-832d-ea1e28bc413d
---
 .../deal.II/include/dofs/dof_constraints.h    |   8 +
 deal.II/deal.II/include/dofs/dof_tools.h      | 113 ++++--
 .../deal.II/source/dofs/dof_constraints.cc    | 256 +++++++++---
 deal.II/deal.II/source/dofs/dof_tools.cc      |  74 +++-
 deal.II/doc/news/2001/c-3-1.html              |  21 +
 .../include/lac/compressed_sparsity_pattern.h | 366 ++++++++++++++++++
 deal.II/lac/include/lac/sparsity_pattern.h    |  10 +
 .../lac/source/compressed_sparsity_pattern.cc | 219 +++++++++++
 deal.II/lac/source/sparsity_pattern.cc        |  24 +-
 tests/deal.II/Makefile.in                     |   3 +-
 tests/deal.II/matrices.cc                     |   4 +-
 tests/deal.II/sparsity_pattern.cc             | 175 +++++++++
 12 files changed, 1164 insertions(+), 109 deletions(-)
 create mode 100644 deal.II/lac/include/lac/compressed_sparsity_pattern.h
 create mode 100644 deal.II/lac/source/compressed_sparsity_pattern.cc
 create mode 100644 tests/deal.II/sparsity_pattern.cc

diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h
index ea0a522b89..9cca92fac6 100644
--- a/deal.II/deal.II/include/dofs/dof_constraints.h
+++ b/deal.II/deal.II/include/dofs/dof_constraints.h
@@ -20,6 +20,7 @@
 #include <base/subscriptor.h>
 
 class SparsityPattern;
+class CompressedSparsityPattern;
 class BlockSparsityPattern;
 template <typename number> class SparseMatrix;
 template <typename number> class BlockSparseMatrix;
@@ -282,6 +283,13 @@ class ConstraintMatrix : public Subscriptor
 				      * patterns.
 				      */
     void condense (BlockSparsityPattern &sparsity) const;
+
+				     /**
+				      * Same function as above, but
+				      * condenses square compressed
+				      * sparsity patterns.
+				      */
+    void condense (CompressedSparsityPattern &sparsity) const;
     
 				     /**
 				      * Condense a given matrix. The associated
diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h
index 7d537ca869..aca12523be 100644
--- a/deal.II/deal.II/include/dofs/dof_tools.h
+++ b/deal.II/deal.II/include/dofs/dof_tools.h
@@ -192,33 +192,44 @@ class DoFTools
 				      * The actual type of the
 				      * sparsity pattern may be
 				      * @ref{SparsityPattern},
-				      * @ref{BlockSparsityPattern}, or any
-				      * other class that satisfies
-				      * similar requirements.
+				      * @ref{CompressedSparsityPattern},
+				      * @ref{BlockSparsityPattern}, or
+				      * any other class that satisfies
+				      * similar requirements. It is
+				      * assumed that the size of the
+				      * sparsity pattern matches the
+				      * number of degrees of freedom
+				      * and that enough unused nonzero
+				      * entries are left to fill the
+				      * sparsity pattern. The nonzero
+				      * entries generated by this
+				      * function are overlaid to
+				      * possible previous content of
+				      * the object, that is previously
+				      * added entries are not deleted.
 				      */
     template <int dim, class SparsityPattern>
-    static void make_sparsity_pattern (const DoFHandler<dim> &dof,
-				       SparsityPattern       &sparsity_pattern);
+    static
+    void
+    make_sparsity_pattern (const DoFHandler<dim> &dof,
+			   SparsityPattern       &sparsity_pattern);
 
 				     /**
 				      * Locate non-zero entries for
-				      * mixed methods.  This function
-				      * does mostly the same as the
-				      * other @p{make_sparsity_pattern},
-				      * but it is specialized for
-				      * mixed finite elements and
-				      * allows to specify which
-				      * variables couple in which
-				      * equation. For example, if
-				      * wanted to solve the Stokes
-				      * equations,
-				      *
-				      *
+				      * ector valued finite elements.
+				      * This function does mostly the
+				      * same as the previous
+				      * @p{make_sparsity_pattern}, but
+				      * it is specialized for vector
+				      * finite elements and allows to
+				      * specify which variables couple
+				      * in which equation. For
+				      * example, if wanted to solve
+				      * the Stokes equations,
 				      * @begin{verbatim}
 				      * -\Delta \vec u + \nabla p = 0,
 				      * \div u                    = 0
 				      * @end{verbatim}
-				      *
 				      * in two space dimensions,
 				      * using stable Q2/Q1 mixed
 				      * elements (using the @ref{FESystem}
@@ -228,7 +239,6 @@ class DoFTools
 				      * rather may want to give the
 				      * following pattern of
 				      * couplings:
-				      *
 				      * @begin{verbatim}
 				      *   1 0 1
 				      *   0 1 1
@@ -263,15 +273,18 @@ class DoFTools
 				      * The actual type of the
 				      * sparsity pattern may be
 				      * @ref{SparsityPattern},
-				      * @ref{BlockSparsityPattern}, or any
-				      * other class that satisfies
+				      * @ref{CompressedSparsityPattern},
+				      * @ref{BlockSparsityPattern}, or
+				      * any other class that satisfies
 				      * similar requirements.
 				      */
-    template<int dim, class SparsityPattern>
-    static void make_sparsity_pattern (const DoFHandler<dim>       &dof,
-				       const std::vector<std::vector<bool> > &mask,
-				       SparsityPattern             &sparsity_pattern);
-
+    template <int dim, class SparsityPattern>
+    static
+    void
+    make_sparsity_pattern (const DoFHandler<dim>                 &dof,
+			   const std::vector<std::vector<bool> > &mask,
+			   SparsityPattern                       &sparsity_pattern);
+    
     				     /**
 				      * Write the sparsity structure
 				      * of the matrix composed of the
@@ -293,12 +306,26 @@ class DoFTools
 				      * matrix yourself, using
 				      * @ref{SparsityPattern}@p{::compress()}.
 				      *
-				      * Since this function is
-				      * obviously useless in one
-				      * spatial dimension, it is not
-				      * implemented.
+				      * The actual type of the
+				      * sparsity pattern may be
+				      * @ref{SparsityPattern},
+				      * @ref{CompressedSparsityPattern},
+				      * @ref{BlockSparsityPattern}, or
+				      * any other class that satisfies
+				      * similar requirements. It is
+				      * assumed that the size of the
+				      * sparsity pattern matches the
+				      * number of degrees of freedom
+				      * and that enough unused nonzero
+				      * entries are left to fill the
+				      * sparsity pattern. The nonzero
+				      * entries generated by this
+				      * function are overlaid to
+				      * possible previous content of
+				      * the object, that is previously
+				      * added entries are not deleted.
 				      */
-    template<int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<dim>           &dof,
 				    const std::vector<unsigned int> &dof_to_boundary_mapping,
@@ -306,8 +333,13 @@ class DoFTools
 
 				     /**
 				      * Declaration of same function
-				      * for different space dimension.
+				      * for one space dimension. The
+				      * resulting sparsity pattern
+				      * couples all dofs of the
+				      * specified end points to each
+				      * other.
 				      */
+    template <class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<1>             &dof,
 				    const std::vector<unsigned int> &dof_to_boundary_mapping,
@@ -340,22 +372,27 @@ class DoFTools
 				      * the function pointers to null
 				      * pointers.
 				      *
-				      * Since this function is
-				      * obviously useless in one
-				      * spatial dimension, it is not
-				      * implemented.
+				      * For the type of the sparsity
+				      * pattern, the same holds as
+				      * said above.
 				      */
-    template<int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<dim> &dof,
 				    const typename FunctionMap<dim>::type &boundary_indicators,
 				    const std::vector<unsigned int>  &dof_to_boundary_mapping,
 				    SparsityPattern    &sparsity); 
 
+
 				     /**
 				      * Declaration of same function
-				      * for different space dimension.
+				      * for one space dimension. The
+				      * resulting sparsity pattern
+				      * couples all dofs of the
+				      * specified end points to each
+				      * other.
 				      */
+    template <class SparsityPattern>
     static void
     make_boundary_sparsity_pattern (const DoFHandler<1>             &dof,
 				    const FunctionMap<1>::type      &boundary_indicators,
diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc
index 9538c3407c..e3cc3f26bd 100644
--- a/deal.II/deal.II/source/dofs/dof_constraints.cc
+++ b/deal.II/deal.II/source/dofs/dof_constraints.cc
@@ -17,6 +17,7 @@
 
 #include <base/memory_consumption.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 #include <lac/vector.h>
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
@@ -423,60 +424,227 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const
   for (unsigned int row=0; row<n_rows; ++row)
     {
       if (distribute[row] == -1)
-					 // regular line. loop over cols
-	for (unsigned int j=sparsity.get_rowstart_indices()[row];
-	     j<sparsity.get_rowstart_indices()[row+1]; ++j)
-	  {
-	    const unsigned int column = sparsity.get_column_numbers()[j];
-	    
+	{
+					   // regular line. loop over
+					   // cols. note that this
+					   // changes the line we are
+					   // presently working on: we
+					   // add additional
+					   // entries. these are put to
+					   // the end of the
+					   // row. however, as
+					   // constrained nodes cannot
+					   // be constrained to other
+					   // constrained nodes, nothing
+					   // will happen if we run into
+					   // these added nodes, as they
+					   // can't be distributed
+					   // further. we might store
+					   // the position of the last
+					   // old entry and stop work
+					   // there, but since operating
+					   // on the newly added ones
+					   // only takes two comparisons
+					   // (column index valid,
+					   // distribute[column]
+					   // necessarily ==-1), it is
+					   // cheaper to not do so and
+					   // run right until the end of
+					   // the line
+	  const unsigned int row_end = sparsity.get_rowstart_indices()[row+1];
+	  for (unsigned int j=sparsity.get_rowstart_indices()[row]; j<row_end; ++j)
+	    {
+	      const unsigned int column = sparsity.get_column_numbers()[j];
+	      
+					       // end of row reached?
+	      if (column == SparsityPattern::invalid_entry)
+		break;
+	      else
+		if (distribute[column] != -1)
+		  {
+						     // distribute entry
+						     // at regular row
+						     // @p{row} and
+						     // irregular column
+						     // sparsity.colnums[j]
+		    for (unsigned int q=0;
+			 q!=lines[distribute[column]].entries.size();
+			 ++q) 
+		      sparsity.add (row,
+				    lines[distribute[column]].entries[q].first);
+		  };
+	    };
+	}
+      else
+					 // row must be
+					 // distributed. note that
+					 // here the present row is
+					 // not touched (unlike above)
+	{
+	  const unsigned int row_end = sparsity.get_rowstart_indices()[row+1];
+	  for (unsigned int j=sparsity.get_rowstart_indices()[row]; j<row_end; ++j)
 					     // end of row reached?
-	    if (column == SparsityPattern::invalid_entry)
+	    if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
 	      break;
 	    else
-	      if (distribute[column] != -1)
-		{
-						   // distribute entry at regular
-						   // row @p{row} and irregular column
+	      {
+		if (distribute[sparsity.get_column_numbers()[j]] == -1)
+						   // distribute entry at irregular
+						   // row @p{row} and regular column
 						   // sparsity.colnums[j]
 		  for (unsigned int q=0;
-		       q!=lines[distribute[column]].entries.size();
-		       ++q) 
-		    sparsity.add (row,
-				  lines[distribute[column]].entries[q].first);
-		};
+		       q!=lines[distribute[row]].entries.size(); ++q) 
+		    sparsity.add (lines[distribute[row]].entries[q].first,
+				  sparsity.get_column_numbers()[j]);
+		else
+						   // distribute entry at irregular
+						   // row @p{row} and irregular column
+						   // sparsity.get_column_numbers()[j]
+		  for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
+		    for (unsigned int q=0;
+			 q!=lines[distribute[sparsity.get_column_numbers()[j]]]
+					.entries.size(); ++q)
+		      sparsity.add (lines[distribute[row]].entries[p].first,
+				    lines[distribute[sparsity.get_column_numbers()[j]]]
+				    .entries[q].first);
+	      };
+	};
+    };
+  
+  sparsity.compress();
+};
+
+
+
+void ConstraintMatrix::condense (CompressedSparsityPattern &sparsity) const
+{
+  Assert (sorted == true, ExcMatrixNotClosed());
+  Assert (sparsity.n_rows() == sparsity.n_cols(),
+	  ExcMatrixNotSquare());
+  
+				   // store for each index whether it
+				   // must be distributed or not. If entry
+				   // is -1, no distribution is necessary.
+				   // otherwise, the number states which
+				   // line in the constraint matrix handles
+				   // this index
+  std::vector<int> distribute(sparsity.n_rows(), -1);
+  
+  for (unsigned int c=0; c<lines.size(); ++c)
+    distribute[lines[c].line] = static_cast<signed int>(c);
+
+  const unsigned int n_rows = sparsity.n_rows();
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+      if (distribute[row] == -1)
+					 // regular line. loop over
+					 // cols. note that as we
+					 // proceed to distribute
+					 // cols, the loop may get
+					 // longer
+	for (unsigned int j=0; j<sparsity.row_length(row); ++j)
+	  {
+	    const unsigned int column = sparsity.column_number(row,j);
+
+ 	    if (distribute[column] != -1)
+	      {
+						 // distribute entry
+						 // at regular row
+						 // @p{row} and
+						 // irregular column
+						 // column. note that
+						 // this changes the
+						 // line we are
+						 // presently working
+						 // on: we add
+						 // additional
+						 // entries. if we add
+						 // another entry at a
+						 // column behind the
+						 // present one, we
+						 // will encounter it
+						 // later on (but
+						 // since it can't be
+						 // further
+						 // constrained, won't
+						 // have to do
+						 // anything about
+						 // it). if we add it
+						 // up front of the
+						 // present column, we
+						 // will find the
+						 // present column
+						 // later on again as
+						 // it was shifted
+						 // back (again
+						 // nothing happens,
+						 // in particular no
+						 // endless loop, as
+						 // when we encounter
+						 // it the second time
+						 // we won't be able
+						 // to add more
+						 // entries as they
+						 // all already exist,
+						 // but we do the same
+						 // work more often
+						 // than necessary,
+						 // and the loop gets
+						 // longer), so move
+						 // the cursor one to
+						 // the right in the
+						 // case that we add
+						 // an entry up front
+						 // that did not exist
+						 // before. check
+						 // whether it existed
+						 // before by tracking
+						 // the length of this
+						 // row
+		unsigned int old_rowlength = sparsity.row_length(row);
+		for (unsigned int q=0;
+		     q!=lines[distribute[column]].entries.size();
+		     ++q) 
+		  {
+		    const unsigned int
+		      new_col = lines[distribute[column]].entries[q].first;
+		    
+		    sparsity.add (row, new_col);
+
+		    const unsigned int new_rowlength = sparsity.row_length(row);
+		    if ((new_col < column) && (old_rowlength != new_rowlength))
+		      ++j;
+		    old_rowlength = new_rowlength;
+		  };
+	      };
 	  }
       else
 					 // row must be distributed
-	for (unsigned int j=sparsity.get_rowstart_indices()[row];
-	     j<sparsity.get_rowstart_indices()[row+1]; ++j)
-					   // end of row reached?
-	  if (sparsity.get_column_numbers()[j] == SparsityPattern::invalid_entry)
-	    break;
-	  else
-	    {
-	      if (distribute[sparsity.get_column_numbers()[j]] == -1)
-						 // distribute entry at irregular
-						 // row @p{row} and regular column
-						 // sparsity.colnums[j]
+	for (unsigned int j=0; j<sparsity.row_length(row); ++j)
+	  {
+	    const unsigned int column = sparsity.column_number(row,j);
+
+	    if (distribute[column] == -1)
+					       // distribute entry at irregular
+					       // row @p{row} and regular column
+					       // sparsity.colnums[j]
+	      for (unsigned int q=0;
+		   q!=lines[distribute[row]].entries.size(); ++q) 
+		sparsity.add (lines[distribute[row]].entries[q].first,
+			      column);
+	    else
+					       // distribute entry at irregular
+					       // row @p{row} and irregular column
+					       // sparsity.get_column_numbers()[j]
+	      for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
 		for (unsigned int q=0;
-		     q!=lines[distribute[row]].entries.size(); ++q) 
-		  sparsity.add (lines[distribute[row]].entries[q].first,
-				sparsity.get_column_numbers()[j]);
-	      else
-						 // distribute entry at irregular
-						 // row @p{row} and irregular column
-						 // sparsity.get_column_numbers()[j]
-		for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
-		  for (unsigned int q=0;
-		       q!=lines[distribute[sparsity.get_column_numbers()[j]]]
-				      .entries.size(); ++q)
-		    sparsity.add (lines[distribute[row]].entries[p].first,
-				  lines[distribute[sparsity.get_column_numbers()[j]]]
-				  .entries[q].first);
-	    };
+		     q!=lines[distribute[sparsity.column_number(row,j)]]
+				    .entries.size(); ++q)
+		  sparsity.add (lines[distribute[row]].entries[p].first,
+				lines[distribute[sparsity.column_number(row,j)]]
+				.entries[q].first);
+	  };
     };
-  
-  sparsity.compress();
 };
 
 
diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc
index 71268fd57a..4ace8c073d 100644
--- a/deal.II/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/deal.II/source/dofs/dof_tools.cc
@@ -27,6 +27,7 @@
 #include <fe/fe_values.h>
 #include <dofs/dof_tools.h>
 #include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 #include <lac/block_sparsity_pattern.h>
 #include <lac/vector.h>
 
@@ -68,6 +69,7 @@ DoFTools::make_sparsity_pattern (const DoFHandler<dim> &dof,
 }
 
 
+
 template <int dim, class SparsityPattern>
 void
 DoFTools::make_sparsity_pattern (const DoFHandler<dim>                 &dof,
@@ -115,8 +117,10 @@ DoFTools::make_sparsity_pattern (const DoFHandler<dim>                 &dof,
 }
 
 
+
 #if deal_II_dimension == 1
 
+template <class SparsityPattern>
 void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>             &dof_handler,
 					  const FunctionMap<1>::type      &function_map,
@@ -160,6 +164,7 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1>             &dof_h
 
 
 
+template <class SparsityPattern>
 void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1> &dof_handler,
 					       const std::vector<unsigned int>  &dof_to_boundary_mapping,
 					       SparsityPattern    &sparsity)
@@ -179,7 +184,7 @@ void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<1> &dof_handler,
 #endif
 
 
-template <int dim>
+template <int dim, class SparsityPattern>
 void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
 					  const std::vector<unsigned int>  &dof_to_boundary_mapping,
@@ -192,9 +197,19 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
 	  ExcDimensionMismatch (sparsity.n_rows(), dof.n_boundary_dofs()));
   Assert (sparsity.n_cols() == dof.n_boundary_dofs(),
 	  ExcDimensionMismatch (sparsity.n_cols(), dof.n_boundary_dofs()));
-  Assert (*max_element(dof_to_boundary_mapping.begin(),
-		       dof_to_boundary_mapping.end()) == sparsity.n_rows()-1,
-	  ExcInternalError());
+#ifdef DEBUG
+  if (true)
+    {
+      unsigned int max_element = 0;
+      for (std::vector<unsigned int>::const_iterator i=dof_to_boundary_mapping.begin();
+	   i!=dof_to_boundary_mapping.end(); ++i)
+	if ((*i != DoFHandler<dim>::invalid_dof_index) &&
+	    (*i > max_element))
+	  max_element = *i;
+      Assert (max_element  == sparsity.n_rows()-1,
+	      ExcInternalError());
+    };
+#endif
 
   const unsigned int dofs_per_face = dof.get_fe().dofs_per_face;
   std::vector<unsigned int> dofs_on_this_face(dofs_per_face);
@@ -216,12 +231,6 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
       {
 	face->get_dof_indices (dofs_on_this_face);
 
-					 // make sure all dof indices have a
-					 // boundary index
-	Assert (*min_element(dofs_on_this_face.begin(),
-			     dofs_on_this_face.end()) != DoFHandler<dim>::invalid_dof_index,
-		ExcInternalError());
-	
 					 // make sparsity pattern for this cell
 	for (unsigned int i=0; i<dofs_per_face; ++i)
 	  for (unsigned int j=0; j<dofs_per_face; ++j) 
@@ -231,7 +240,8 @@ DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
 };
 
 
-template <int dim>
+
+template <int dim, class SparsityPattern>
 void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
 					       const typename FunctionMap<dim>::type  &boundary_indicators,
 					       const std::vector<unsigned int>  &dof_to_boundary_mapping,
@@ -270,11 +280,6 @@ void DoFTools::make_boundary_sparsity_pattern (const DoFHandler<dim>& dof,
       {
 	face->get_dof_indices (dofs_on_this_face);
 
-					 // make sure all dof indices have a
-					 // boundary index
-	Assert (*min_element(dofs_on_this_face.begin(),
-			     dofs_on_this_face.end()) != DoFHandler<dim>::invalid_dof_index,
-		ExcInternalError());
 					 // make sparsity pattern for this cell
 	for (unsigned int i=0; i<dofs_per_face; ++i)
 	  for (unsigned int j=0; j<dofs_per_face; ++j)
@@ -2163,35 +2168,64 @@ DoFTools::make_flux_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
 				      SparsityPattern    &,
 				      const FullMatrix<double>&,
 				      const FullMatrix<double>&);
+#endif
+
 template void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
 					  const std::vector<unsigned int>  &,
 					  SparsityPattern    &);
 template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+					  const std::vector<unsigned int>  &,
+					  CompressedSparsityPattern    &);
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+					  const std::vector<unsigned int>  &,
+					  BlockSparsityPattern    &);
+template void
 DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
 					  const FunctionMap<deal_II_dimension>::type  &boundary_indicators,
 					  const std::vector<unsigned int>  &dof_to_boundary_mapping,
 					  SparsityPattern    &sparsity);
-#endif
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+					  const FunctionMap<deal_II_dimension>::type  &boundary_indicators,
+					  const std::vector<unsigned int>  &dof_to_boundary_mapping,
+					  CompressedSparsityPattern    &sparsity);
+template void
+DoFTools::make_boundary_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+					  const FunctionMap<deal_II_dimension>::type  &boundary_indicators,
+					  const std::vector<unsigned int>  &dof_to_boundary_mapping,
+					  BlockSparsityPattern    &sparsity);
 
 template void
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
 				 SparsityPattern    &sparsity);
 
+template void
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
+				 CompressedSparsityPattern    &sparsity);
+
 template void
 DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
 				 BlockSparsityPattern                &sparsity);
 
 template void 
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
 				 const std::vector<std::vector<bool> > &mask,
 				 SparsityPattern             &sparsity);
 
 template void 
-DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension>& dof,
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
+				 const std::vector<std::vector<bool> > &mask,
+				 CompressedSparsityPattern             &sparsity);
+
+template void 
+DoFTools::make_sparsity_pattern (const DoFHandler<deal_II_dimension> &dof,
 				 const std::vector<std::vector<bool> > &mask,
 				 BlockSparsityPattern        &sparsity);
 
+
 template
 void
 DoFTools::distribute_cell_to_dof_vector (const DoFHandler<deal_II_dimension> &dof_handler,
diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index 4f3d503521..2a35c76891 100644
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -305,6 +305,16 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       New: Class <code class="class">CompressedSparsityPattern</code>
+       may be used as an intermediate form of the <code
+       class="class">SparsityPattern</code> class if memory
+       requirements are tight during construction of the sparsity
+       pattern.
+       <br>
+       (WB 2001/06/22)
+       </p>
+
   <li> <p>
        New: There are now functions 
        <code class="member">SparsityPattern::copy_from</code> and
@@ -441,6 +451,17 @@ documentation, etc</a>.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p>
+       New: Implement
+       <code class="member">DoFTools::make_sparsity_pattern</class>,
+       <code class="member">DoFTools::make_boundary_sparsity_pattern</class>,
+       and 
+       <code class="member">ConstraintMatrix::condense</class> to work on
+       the <code class="class">CompressedSparsityPattern</code> class.
+       <br>
+       (WB 2001/06/22)
+       </p>
+
   <li> <p>
        New: <code class="member">FE*Values::get_quadrature</class>
        returns a reference to the quadrature formula used by a
diff --git a/deal.II/lac/include/lac/compressed_sparsity_pattern.h b/deal.II/lac/include/lac/compressed_sparsity_pattern.h
new file mode 100644
index 0000000000..c5212885e5
--- /dev/null
+++ b/deal.II/lac/include/lac/compressed_sparsity_pattern.h
@@ -0,0 +1,366 @@
+//----------------------------  compressed_sparsity_pattern.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  compressed_sparsity_pattern.h  ---------------------------
+#ifndef __deal2__compressed_sparsity_pattern_h
+#define __deal2__compressed_sparsity_pattern_h
+
+
+#include <base/exceptions.h>
+#include <base/subscriptor.h>
+
+template <typename number> class SparseMatrix;
+
+#include <vector>
+#include <set>
+
+
+/**
+ * This class acts as an intermediate form of the
+ * @ref{SparsityPattern} class. From the interface it mostly
+ * represents a @ref{SparsityPattern} object that is kept compressed
+ * at all times. However, since the final sparsity pattern is not
+ * known while constructing it, keeping the pattern compressed at all
+ * times can only be achieved at the expense of either increased
+ * memory or run time consumption upon use. The main purpose of this
+ * class is to avoid some memory bottlenecks, so we chose to implement
+ * it memory conservative, but the chosen data format is too unsuited
+ * to be used for actual matrices. It is therefore necessary to first
+ * copy the data of this object over to an object of type
+ * @ref{SparsityPattern} before using it in actual matrices.
+ *
+ * Another viewpoint is that this class does not need up front
+ * allocation of a certain amount of memory, but grows as necessary.
+ *
+ *
+ * @sect3{Rationale}
+ *
+ * When constructing the sparsity pattern of a matrix, you usually
+ * first have to provide an empty sparsity pattern object with a fixed
+ * maximal number of entries per row. To find out about this maximal
+ * row length, one usually calls the function
+ * @ref{DoFHandler}@p{::max_couplings_per_dof} which returns an
+ * estimate for that quantity. While this estimate is usually quite
+ * good in 2d and exact in 1d, it is often significantly too large in
+ * 3d and especially for higher order elements. Furthermore, normally
+ * only a small fraction of the rows of a matrix will end up having
+ * the maximal number of nonzero entries per row (usually those nodes
+ * adjacent to hanging nodes), most have much less. In effect, the
+ * empty @ref{SparsityPattern} object has allocated much too much
+ * memory. Although this unnecessarily allocated memory is later freed
+ * when @ref{SparsityPattern}@p{::compress} is called, this
+ * overallocation has, with higher order elements and in 3d, sometimes
+ * been so large that the program aborted due to lack of memory.
+ *
+ * This class therefore provides an alternative representation of a
+ * sparsity pattern: we don't specify a maximal row length initially,
+ * but store a set of column indices indicating possible nonzero
+ * entries in the sparsity pattern for each row. This is very much
+ * like the final "compressed" format used in the
+ * @ref{SparsityPattern} object after compression, but uses a less
+ * compact memory storage format, since the exact number of entries
+ * per row is only known a posteriori and since it may change (for the
+ * @ref{SparsityPattern} class, no more changes are allowed after
+ * compressing it). We can therefore not store all the column indices
+ * in a big array, but have to use a vector of sets. This can later be
+ * used to actually initialize a @ref{SparsityPattern} object with the
+ * then final set of necessary indices.
+ *
+ *
+ * @sect3{Interface}
+ *
+ * Since this class is intended as an intermediate replacement of the
+ * @ref{SparsityPattern} class, it has mostly the same interface, with
+ * small changes where necessary. In particular, the @ref{add}
+ * function, and the functions inquiring properties of the sparsity
+ * pattern are the same.
+ *
+ *
+ * @sect3{Usage}
+ *
+ * Use this class as follows:
+ * @begin{verbatim}
+ * CompressedSparsityPattern compressed_pattern (dof_handler.n_dofs());
+ * DoFTools::make_sparsity_pattern (dof_handler,
+ *                                  compressed_pattern);
+ * constraints.condense (compressed_pattern);
+ *
+ * SparsityPattern sp;
+ * sp.copy_from (compressed_pattern);
+ * @end{verbatim}
+ *
+ *
+ * @author Wolfgang Bangerth, 2001
+ */
+class CompressedSparsityPattern : public Subscriptor
+{
+  public:
+				     /**
+				      * Initialize the matrix empty,
+				      * that is with no memory
+				      * allocated. This is useful if
+				      * you want such objects as
+				      * member variables in other
+				      * classes. You can make the
+				      * structure usable by calling
+				      * the @p{reinit} function.
+				      */
+    CompressedSparsityPattern ();
+    
+				     /**
+				      * Copy constructor. This constructor is
+				      * only allowed to be called if the matrix
+				      * structure to be copied is empty. This is
+				      * so in order to prevent involuntary
+				      * copies of objects for temporaries, which
+				      * can use large amounts of computing time.
+				      * However, copy constructors are needed
+				      * if yo want to use the STL data types
+				      * on classes like this, e.g. to write
+				      * such statements like
+				      * @p{v.push_back (CompressedSparsityPattern());},
+				      * with @p{v} a vector of @p{CompressedSparsityPattern}
+				      * objects.
+				      *
+				      * Usually, it is sufficient to
+				      * use the explicit keyword to
+				      * disallow unwanted temporaries,
+				      * but for the STL vectors, this
+				      * does not work. Since copying a
+				      * structure like this is not
+				      * useful anyway because multiple
+				      * matrices can use the same
+				      * sparsity structure, copies are
+				      * only allowed for empty
+				      * objects, as described above.
+				      */
+    CompressedSparsityPattern (const CompressedSparsityPattern &);
+
+				     /**
+				      * Initialize a rectangular
+				      * matrix with @p{m} rows and
+				      * @p{n} columns.
+				      */
+    CompressedSparsityPattern (const unsigned int m,
+			       const unsigned int n);
+    
+				     /**
+				      * Initialize a square matrix of
+				      * dimension @p{n}.
+				      */
+    CompressedSparsityPattern (const unsigned int n);
+
+				     /**
+				      * Copy operator. For this the
+				      * same holds as for the copy
+				      * constructor: it is declared,
+				      * defined and fine to be called,
+				      * but the latter only for empty
+				      * objects.
+				      */
+    CompressedSparsityPattern & operator = (const CompressedSparsityPattern &);
+    
+				     /**
+				      * Reallocate memory and set up
+				      * data structures for a new
+				      * matrix with @p{m} rows and
+				      * @p{n} columns, with at most
+				      * @p{max_per_row} nonzero
+				      * entries per row.
+				      */
+    void reinit (const unsigned int m,
+		 const unsigned int n);
+    
+				     /**
+				      * Since this object is kept
+				      * compressed at all times anway,
+				      * this function does nothing,
+				      * but is declared to make the
+				      * interface of this class as
+				      * much alike as that of the
+				      * @ref{SparsityPattern} class.
+				      */
+    void compress ();
+    
+				     /**
+				      * Return whether the object is
+				      * empty. It is empty if no
+				      * memory is allocated, which is
+				      * the same as that both
+				      * dimensions are zero.
+				      */
+    bool empty () const;
+
+				     /**
+				      * Return the maximum number of
+				      * entries per row. Note that
+				      * this number may change as
+				      * entries are added.
+				      */
+    unsigned int max_entries_per_row () const;
+
+				     /**
+				      * Add a nonzero entry to the
+				      * matrix. If the entry already
+				      * exists, nothing bad happens.
+				      */
+    void add (const unsigned int i, 
+	      const unsigned int j);
+    
+				     /**
+				      * Make the sparsity pattern
+				      * symmetric by adding the
+				      * sparsity pattern of the
+				      * transpose object.
+				      *
+				      * This function throws an
+				      * exception if the sparsity
+				      * pattern does not represent a
+				      * square matrix.
+				      */
+    void symmetrize ();
+    
+				     /**
+				      * Print the sparsity of the matrix
+				      * in a format that @p{gnuplot} understands
+				      * and which can be used to plot the
+				      * sparsity pattern in a graphical
+				      * way. The format consists of pairs
+				      * @p{i j} of nonzero elements, each
+				      * representing one entry of this
+				      * matrix, one per line of the output
+				      * file. Indices are counted from
+				      * zero on, as usual. Since sparsity
+				      * patterns are printed in the same
+				      * way as matrices are displayed, we
+				      * print the negative of the column
+				      * index, which means that the
+				      * @p{(0,0)} element is in the top left
+				      * rather than in the bottom left
+				      * corner.
+				      *
+				      * Print the sparsity pattern in
+				      * gnuplot by setting the data style
+				      * to dots or points and use the
+				      * @p{plot} command.
+				      */
+    void print_gnuplot (std::ostream &out) const;
+
+				     /**
+				      * Return number of rows of this
+				      * matrix, which equals the dimension
+				      * of the image space.
+				      */
+    unsigned int n_rows () const;
+
+				     /**
+				      * Return number of columns of this
+				      * matrix, which equals the dimension
+				      * of the range space.
+				      */
+    unsigned int n_cols () const;
+
+				     /**
+				      * Number of entries in a specific row.
+				      */
+    unsigned int row_length (const unsigned int row) const;
+
+				     /**
+				      * Access to column number field.
+				      * Return the column number of
+				      * the @p{index}th entry in @p{row}.
+				      */
+    unsigned int column_number (const unsigned int row,
+				const unsigned int index) const;
+
+				     /**
+				      * Compute the bandwidth of the matrix
+				      * represented by this structure. The
+				      * bandwidth is the maximum of
+				      * $|i-j|$ for which the index pair
+				      * $(i,j)$ represents a nonzero entry
+				      * of the matrix.
+				      */
+    unsigned int bandwidth () const;
+
+				     /**
+				      * Return the number of nonzero elements of
+				      * this matrix. Actually, it returns the
+				      * number of entries in the sparsity
+				      * pattern; if any of the entries should
+				      * happen to be zero, it is counted
+				      * anyway.
+				      *
+				      * This function may only be called if the
+				      * matrix struct is compressed. It does not
+				      * make too much sense otherwise anyway.
+				      */
+    unsigned int n_nonzero_elements () const;
+
+    				     /**
+				      * Exception
+				      */
+    DeclException2 (ExcInvalidIndex,
+		    int, int,
+		    << "The given index " << arg1
+		    << " should be less than " << arg2 << ".");
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcInvalidConstructorCall);
+				     /**
+				      * Exception
+				      */
+    DeclException0 (ExcNotSquare);
+    
+  private:
+				     /**
+				      * Number of rows that this sparsity
+				      * structure shall represent.
+				      */
+    unsigned int rows;
+
+				     /**
+				      * Number of columns that this sparsity
+				      * structure shall represent.
+				      */
+    unsigned int cols;
+
+				     /**
+				      * Actual data: store for each
+				      * row the set of nonzero
+				      * entries.
+				      */
+    std::vector<std::set<unsigned int> > column_indices;
+
+    friend class SparsityPattern;
+};
+
+
+/*---------------------- Inline functions -----------------------------------*/
+
+
+inline
+unsigned int
+CompressedSparsityPattern::n_rows () const
+{
+  return rows;
+};
+
+
+inline
+unsigned int
+CompressedSparsityPattern::n_cols () const
+{
+  return cols;
+};
+
+#endif
diff --git a/deal.II/lac/include/lac/sparsity_pattern.h b/deal.II/lac/include/lac/sparsity_pattern.h
index d4baa0cdf3..2940dbf3fc 100644
--- a/deal.II/lac/include/lac/sparsity_pattern.h
+++ b/deal.II/lac/include/lac/sparsity_pattern.h
@@ -18,6 +18,7 @@
 #include <base/subscriptor.h>
 
 template <typename number> class SparseMatrix;
+class CompressedSparsityPattern;
 
 #include <vector>
 #include <iterator>
@@ -460,6 +461,15 @@ class SparsityPattern : public Subscriptor
 		    const unsigned int    n_cols,
 		    const ForwardIterator begin,
 		    const ForwardIterator end);
+
+				     /**
+				      * Copy data from an object of
+				      * type
+				      * @ref{CompressedSparsityPattern}. Previous
+				      * content of this object is
+				      * lost.
+				      */
+    void copy_from (const CompressedSparsityPattern &csp);
     
 				     /**
 				      * Return whether the object is empty. It
diff --git a/deal.II/lac/source/compressed_sparsity_pattern.cc b/deal.II/lac/source/compressed_sparsity_pattern.cc
new file mode 100644
index 0000000000..34f27f8f8f
--- /dev/null
+++ b/deal.II/lac/source/compressed_sparsity_pattern.cc
@@ -0,0 +1,219 @@
+//----------------------------  compressed_sparsity_pattern.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  compressed_sparsity_pattern.cc  ---------------------------
+
+
+#include <lac/compressed_sparsity_pattern.h>
+
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+#include <functional>
+
+CompressedSparsityPattern::CompressedSparsityPattern () :
+		rows(0),
+		cols(0)
+{};
+
+
+
+CompressedSparsityPattern::CompressedSparsityPattern (const CompressedSparsityPattern &s) :
+		Subscriptor(),
+		rows(0),
+		cols(0)
+{
+  Assert (s.rows == 0, ExcInvalidConstructorCall());
+  Assert (s.cols == 0, ExcInvalidConstructorCall());
+};
+
+
+
+CompressedSparsityPattern::CompressedSparsityPattern (const unsigned int m,
+						      const unsigned int n) 
+		: rows(0),
+		  cols(0)
+{
+  reinit (m,n);
+};
+
+
+
+CompressedSparsityPattern::CompressedSparsityPattern (const unsigned int n)
+		: rows(0),
+		  cols(0)
+{
+  reinit (n,n);
+};
+
+
+
+CompressedSparsityPattern &
+CompressedSparsityPattern::operator = (const CompressedSparsityPattern &s)
+{
+  Assert (s.rows == 0, ExcInvalidConstructorCall());
+  Assert (s.cols == 0, ExcInvalidConstructorCall());
+
+  Assert (rows == 0, ExcInvalidConstructorCall());
+  Assert (cols == 0, ExcInvalidConstructorCall());
+
+  return *this;
+};
+
+
+
+void
+CompressedSparsityPattern::reinit (const unsigned int m,
+				   const unsigned int n)
+{
+  rows = m;
+  cols = n;
+
+  std::vector<std::set<unsigned int> > new_column_indices (rows);
+  column_indices.swap (new_column_indices);
+}
+
+
+
+void
+CompressedSparsityPattern::compress ()
+{};
+
+
+
+bool
+CompressedSparsityPattern::empty () const
+{
+  return ((rows==0) && (cols==0));
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::max_entries_per_row () const 
+{
+  unsigned int m = 0;
+  for (unsigned int i=1; i<rows; ++i)
+    m = std::max (m, column_indices[i].size());
+
+  return m;
+};
+
+
+
+void
+CompressedSparsityPattern::add (const unsigned int i,
+				const unsigned int j)
+{
+  Assert (i<rows, ExcInvalidIndex(i,rows));
+  Assert (j<cols, ExcInvalidIndex(j,cols));
+
+				   // the std::set automatically
+				   // assures uniqueness and
+				   // sortedness of the column indices
+  column_indices[i].insert (j);
+};
+
+
+
+void
+CompressedSparsityPattern::symmetrize () 
+{
+  Assert (rows==cols, ExcNotSquare());
+
+				   // loop over all elements presently
+				   // in the sparsity pattern and add
+				   // the transpose element. note:
+				   //
+				   // 1. that the sparsity pattern
+				   // changes which we work on, but
+				   // not the present row
+				   //
+				   // 2. that the @p{add} function can
+				   // be called on elements that
+				   // already exist without any harm
+  for (unsigned int row=0; row<rows; ++row)
+    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
+	 i!=column_indices[row].end(); ++i)
+				       // add the transpose entry if
+				       // this is not the diagonal
+      if (row != *i)
+	add (*i, row);
+};
+
+
+
+void
+CompressedSparsityPattern::print_gnuplot (std::ostream &out) const
+{
+  for (unsigned int row=0; row<rows; ++row)
+    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
+	 i!=column_indices[row].end(); ++i)
+				       // while matrix entries are
+				       // usually written (i,j),
+				       // with i vertical and j
+				       // horizontal, gnuplot output
+				       // is x-y, that is we have to
+				       // exchange the order of
+				       // output
+      out << *i << " " << -static_cast<signed int>(row) << std::endl;
+
+  AssertThrow (out, ExcIO());
+}
+
+
+
+unsigned int
+CompressedSparsityPattern::row_length (const unsigned int row) const
+{
+  return column_indices[row].size();
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::column_number (const unsigned int row,
+					  const unsigned int index) const
+{
+  Assert (index < column_indices[row].size(),
+	  ExcIndexRange (index, 0, column_indices[row].size()));
+  std::set<unsigned int>::const_iterator p = column_indices[row].begin();
+  std::advance (p, index);
+  return *p;
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::bandwidth () const
+{
+  unsigned int b=0;
+  for (unsigned int row=0; row<rows; ++row)
+    for (std::set<unsigned int>::const_iterator i=column_indices[row].begin();
+	 i!=column_indices[row].end(); ++i)
+      if (static_cast<unsigned int>(abs(static_cast<int>(row-*i))) > b)
+	b = abs(static_cast<signed int>(row-*i));
+
+  return b;
+};
+
+
+
+unsigned int
+CompressedSparsityPattern::n_nonzero_elements () const
+{
+  unsigned int n=0;
+  for (unsigned int i=0; i<rows; ++i)
+    n += column_indices[i].size();
+  return n;
+};
diff --git a/deal.II/lac/source/sparsity_pattern.cc b/deal.II/lac/source/sparsity_pattern.cc
index d7fba38f51..6a1d1512d2 100644
--- a/deal.II/lac/source/sparsity_pattern.cc
+++ b/deal.II/lac/source/sparsity_pattern.cc
@@ -12,7 +12,8 @@
 //----------------------------  sparsity_pattern.cc  ---------------------------
 
 
-#include <lac/sparse_matrix.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
 
 #include <iostream>
 #include <iomanip>
@@ -447,6 +448,16 @@ SparsityPattern::compress ()
 
 
 
+void
+SparsityPattern::copy_from (const CompressedSparsityPattern &csp) 
+{
+  copy_from (csp.n_rows(), csp.n_cols(),
+	     csp.column_indices.begin(),
+	     csp.column_indices.end());
+};
+
+
+
 bool
 SparsityPattern::empty () const
 {
@@ -583,7 +594,8 @@ SparsityPattern::symmetrize ()
 				   // the transpose element. note:
 				   //
 				   // 1. that the sparsity pattern
-				   // changes which we work on
+				   // changes which we work on, but
+				   // not the present row
 				   //
 				   // 2. that the @p{add} function can
 				   // be called on elements that
@@ -599,8 +611,12 @@ SparsityPattern::symmetrize ()
 	  break;
 
 					 // otherwise add the
-					 // transpose entry
-	add (colnums[k], row);
+					 // transpose entry if this is
+					 // not the diagonal (that
+					 // would not harm, only take
+					 // time to check up)
+	if (colnums[k] != row)
+	  add (colnums[k], row);
       };
 };
 
diff --git a/tests/deal.II/Makefile.in b/tests/deal.II/Makefile.in
index 761b1be944..4cb007ae05 100644
--- a/tests/deal.II/Makefile.in
+++ b/tests/deal.II/Makefile.in
@@ -42,12 +42,13 @@ wave-test-3.exe              : wave-test-3.go                        $(lib-2d)
 support_point_map.exe        : support_point_map.go        $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 filtered_matrix.exe          : filtered_matrix.go          $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 boundaries.exe               : boundaries.go               $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
+sparsity_pattern.exe         : sparsity_pattern.go         $(lib-1d) $(lib-2d) $(lib-1d) $(libraries)
 
 tests = grid_test grid_transform dof_test data_out derivatives gradients constraints mg \
         mglocal block_matrices second_derivatives derivative_approximation \
         matrices error_estimator intergrid_constraints intergrid_map \
         wave-test-3 dof_renumbering support_point_map filtered_matrix \
-	boundaries
+	boundaries sparsity_pattern
 
 ############################################################
 
diff --git a/tests/deal.II/matrices.cc b/tests/deal.II/matrices.cc
index 19d526c181..97e2551f2e 100644
--- a/tests/deal.II/matrices.cc
+++ b/tests/deal.II/matrices.cc
@@ -1,4 +1,4 @@
-//----------------------------  derivative_approximation.cc  ---------------------------
+//----------------------------  matrices.cc  ---------------------------
 //    $Id$
 //    Version: $Name$
 //
@@ -9,7 +9,7 @@
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  derivative_approximation.cc  ---------------------------
+//----------------------------  matrices.cc  ---------------------------
 
 
 /* Author: Wolfgang Bangerth, University of Heidelberg, 2001 */
diff --git a/tests/deal.II/sparsity_pattern.cc b/tests/deal.II/sparsity_pattern.cc
new file mode 100644
index 0000000000..ec18708351
--- /dev/null
+++ b/tests/deal.II/sparsity_pattern.cc
@@ -0,0 +1,175 @@
+//----------------------------  sparsity_pattern.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2000, 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  sparsity_pattern.cc  ---------------------------
+
+
+/* Author: Wolfgang Bangerth, University of Heidelberg, 2001 */
+
+// check that the direct generation of the sparsity pattern and that
+// via the CompressedSparsityPattern result in the same
+
+
+#include <base/logstream.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/fe_q.h>
+#include <fe/fe_system.h>
+
+#include <fstream>
+
+
+
+bool operator == (const SparsityPattern &sp1,
+		  const SparsityPattern &sp2)
+{
+  if (sp1.n_nonzero_elements() != sp2.n_nonzero_elements())
+    return false;
+  
+  for (unsigned int i=0; i<sp1.n_nonzero_elements(); ++i)
+    if (sp1.get_column_numbers()[i] !=
+	sp2.get_column_numbers()[i])
+      return false;
+
+  for (unsigned int i=0; i<sp1.n_rows(); ++i)
+    if (sp1.get_rowstart_indices()[i] !=
+	sp2.get_rowstart_indices()[i])
+      return false;
+  
+  return true;
+};
+
+
+
+template <int dim>
+void
+check_boundary (const DoFHandler<dim> &dof)
+{
+  std::vector<unsigned int> dof_to_boundary_mapping;
+  DoFTools::map_dof_to_boundary_indices (dof,
+					 dof_to_boundary_mapping);
+
+				   // first way: direct generation
+  SparsityPattern sparsity_1(dof.n_boundary_dofs(),
+			     dof.max_couplings_between_boundary_dofs());
+  DoFTools::make_boundary_sparsity_pattern (dof,
+					    dof_to_boundary_mapping,
+					    sparsity_1);
+  sparsity_1.compress ();
+
+				   // second way: via a CompressedSparsityPattern
+  SparsityPattern sparsity_2;
+  CompressedSparsityPattern csp(dof.n_boundary_dofs());
+  DoFTools::make_boundary_sparsity_pattern (dof,
+					    dof_to_boundary_mapping,
+					    csp);
+  sparsity_2.copy_from (csp);
+
+				   // the exact content of sparsity
+				   // patterns is checked in other
+				   // tests, so only make sure that
+				   // sparsity_[12] are equal
+  deallog << __PRETTY_FUNCTION__
+	  << " -- "
+	  << (sparsity_1 == sparsity_2 ? "ok" : "failed")
+	  << std::endl;
+};
+
+
+
+
+template <int dim>
+void
+check ()
+{
+  Triangulation<dim> tr;  
+  if (dim==2)
+    GridGenerator::hyper_ball(tr, Point<dim>(), 1);
+  else
+    GridGenerator::hyper_cube(tr, -1,1);
+  tr.refine_global (1);
+  tr.begin_active()->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  tr.begin_active(2)->set_refine_flag ();
+  tr.execute_coarsening_and_refinement ();
+  if (dim==1)
+    tr.refine_global(2);
+
+				   // create a system element composed
+				   // of one Q1 and one Q2 element
+  FESystem<dim> element(FE_Q<dim>(1), 1,
+			FE_Q<dim>(2), 1);
+  DoFHandler<dim> dof(tr);
+  dof.distribute_dofs(element);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  constraints.close ();
+
+				   // create sparsity pattern. note
+				   // that different components should
+				   // not couple, so use pattern
+  std::vector<std::vector<bool> > mask (2, std::vector<bool>(2, false));
+  mask[0][0] = mask[1][1] = true;
+
+				   // first way: directly
+  SparsityPattern sparsity_1 (dof.n_dofs(), dof.n_dofs());
+  DoFTools::make_sparsity_pattern (dof, mask, sparsity_1);
+  constraints.condense (sparsity_1);
+  sparsity_1.compress ();
+
+				   // second way: via CompressedSparsityPattern
+  SparsityPattern sparsity_2;
+  CompressedSparsityPattern csp (dof.n_dofs());
+  DoFTools::make_sparsity_pattern (dof, mask, csp);
+  constraints.condense (csp);
+  sparsity_2.copy_from (csp);
+
+
+				   // the exact content of sparsity
+				   // patterns is checked in other
+				   // tests, so only make sure that
+				   // sparsity_[12] are equal
+  deallog << __PRETTY_FUNCTION__
+	  << " -- "
+	  << (sparsity_1 == sparsity_2 ? "ok" : "failed")
+	  << std::endl;
+  
+  check_boundary (dof);
+};
+
+
+
+int main ()
+{
+  std::ofstream logfile ("sparsity_pattern.output");
+  logfile.precision (2);
+  logfile.setf(std::ios::fixed);  
+  deallog.attach(logfile);
+  deallog.depth_console (0);
+
+  deallog.push ("1d");
+  check<1> ();
+  deallog.pop ();
+  deallog.push ("2d");
+  check<2> ();
+  deallog.pop ();
+  deallog.push ("3d");
+  check<3> ();
+  deallog.pop ();
+}
-- 
2.39.5